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Abstract 


The  Finite  Element  Method  was  used  to  solve  the  nonlinear  electron 
plasma  equations  for  the  System-Generated  Electromagnetic  Pulse  boundary 
layer  in  one  spatial  dimension.  These  equations  were  solved  in  dis¬ 
tance-velocity  phase  space  using  a  rectangular  finite  element  mesh. 
Linear  approximations  were  used  for  both  the  trial  and  weight  functions 
for  each  element.  The  advection  terms  in  the  Vlasov  plasma  equation 
were  treated  with  the  Heinrich  upwinding  technique. 

The  time  integration  was  performed  using  an  explicit  two-step 
Lax-Wendroff  procedure.  The  system  of  algebraic  equations  was  solved 
with  a  fully-packed  Gauss-Seidel  iteration  scheme. 

A  value  of  2/3  for  the  upwinding  parameter  was  found  to  provide  the 
best  compromise  between  dispersion  of  the  pulse  and  computer  storage 
requirements.  The  savings  in  computer  memory  results  in  increased 
execution  speed  for  the  algorithm.  Also,  it  is  shown  that  the  numerical 
scheme  does  not  permit  spurious  pulse  reflections  from  the  edges  of  the 
mesh. 

Results  for  several  test  cases  are  presented.  Comparisons  are 
given  which  show  favorable  agreement  for  the  finite  element  technique 
with  other  solution  methods. 

Empirical  relationships  for  the  mesh  parameters  are  given  which 
must  be  followed  in  order  to  produce  valid  results  with  the  numerical 
scheme  developed. 
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THE  FINITE  ELEMENT  METHOD  APPLIED 


TO  THE  SYSTEM-GENERATED  ELECTROMAGNETIC  PULSE 
ROUNDARY  LAYER 

I.  Introduction 

DEFINITION  OF  SCEMP.  System-Generated  Electromagnetic  Pulse 
(SGEMP)  is  an  effect  of  a  nuclear  weapon  detonation.  It  will  occur 
whenever  the  ionizing  radiations  (X-rays,  gamma  rays,  and  neutrons)  from 
a  nuclear  event  are  incident  upon  an  object  in  a  low  pressure  (high 
altitude)  environment.  The  radiation  interacts  with  the  object  via  the 
Compton  and  photoelectric  processes.  These  mechanisms  produce  free 
electrons  inside  the  object  and  in  the  region  surrounding  it.  The  elec¬ 
tromagnetic  field  created  by  the  electrons  in  motion  about  the  object  is 
called  SGEMP. 

A  nuclear  burst  in  space  is  of  particular  interest  to  the  Air 
Force.  Line-of-s ight  radiation  from  such  a  burst  will  be  received  by  a 
satellite  unattenuated  by  atmospheric  interactions.  The  line-of-sight 
radiation  arrives  at  the  satellite  affected  only  by  the  spherical  di¬ 
vergence  factor  of  1/4-rtr2,  where  r  =  distance  from  the  burst  point  to 
the  satellite. 

The  depth  of  interaction  into  the  material  of  the  satellite  depends 
on  the  energy  and  type  of  radiation.  High  energy  X-rays,  gamma  rays  and 
neutrons  will  penetrate  a  few,  to  tens  of  centimeters,  into  the  material 
before  interactions  take  place.  These  interactions  lead  to  direct 
injection  of  current  into  the  satellite's  circuitry. 

However,  the  penetrating  radiations  are  only  a  minor  fraction  of  a 
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nuclear  weapon's  energy  output.  Approximately  80%  of  the  energy 
released  from  a  nuclear  hurst  is  in  the  form  of  X-rays  (Pef  1,  Chapter 
7).  A  fraction  of  these  X-rays  are  in  the  energy  range  of  1  keV.  These 
relatively  low-energy  photons  will  primarily  interact  with  the  surface 
atoms  of  the  satellite  (usually  a  metal  such  as  Aluminum).  About  half 
of  all  the  electrons  generated  by  the  photoelectric  effect  (the  major 
interaction  mechanism  at  these  photon  energies)  are  back-scattered  into 
space  off  of  the  satellite  (Ref  2,  Chapters  23-25).  Consequently,  elec¬ 
tromagnetic  fields  will  surround  the  satellite,  generating  surface 
currents  in  the  process.  These  surface  currents  may  also  find  their  way 
into  the  circuitry  of  the  satellite,  creating  a  potential  hazard  to 
electronic  components.  The  creation  of  surface  currents  via  this  pro¬ 
cess  is  known  as  the  external  (or  outside)  System-Generated  Electro¬ 
magnetic  Pulse  problem,  or  simply  SGEMP,  and  is  the  phenomenon  analyzed 
here. 

IMPORTANCE  OF  UNDERSTANDING  SGEMP.  The  ability  of  a  satellite  to 
withstand  a  nuclear  environment  is  critical  to  our  national  defense. 
Therefore,  a  quantitative  knowledge  of  the  nuclear  threat  to  these 
systems  is  necessary.  The  SGEMP  effect  is  an  important  consideration 
for  space  systems.  Realistic  testing  is  one  of  the  best  methods  to 
measure  the  impact  of  SGEMP  on  space  vehicles.  Unfortunately,  the 
creation  of  the  right  kind  of  radiation  with  the  proper  energy,  flux, 
and  time  history  is  difficult  to  produce  in  the  laboratory.  Underground 
nuclear  testing  can  be  done,  but  the  cost,  and  large  size  of  some 
systems  makes  such  experiments  difficult  to  accomplish,  at  best. 

As  a  result,  the  Air  Force  has  instituted  a  technical  program  for 
the  theoretical  understanding  of  SGEMP.  This  program  encompasses  both 
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analytical  and  computational  methods.  Researchers  have  studied  several 
important  properties  of  the  SGEMP  process  since  the  institution  of  this 
program.  One  of  these  is  the  boundary  layer. 

THE  BOUNDARY  LAYER.  The  electron  boundary  layer  is  created  near 
the  surface  of  the  satellite.  This  effect  is  due  to  the  large,  normal 
electric  field  created  by  the  removal  of  charge  from  the  material  of  the 
satellite.  Since  this  electric  field  is  directed  away  from  the  surface, 
electrons  experience  a  force  directed  back  towards  the  satellite.  In 
fact,  the  electric  field  can  get  so  strong  that  subsequently  ejected 
electrons  of  lower  energy  will  not  be  able  to  penetrate  this  small 
region  near  the  satellite.  This  region  will  contain  a  high  density  of 
electrons,  some  which  have  been  slowed  down  considerably,  others  which 
have  actually  turned  around.  It  is  even  possible  for  electrons  to  exist 
in  this  area  under  quasi-static  equilibrium  conditions. 

This  region  of  dense  electrons  which  retards  the  progress  of  other 
emitted  electrons  is  known  as  the  boundary  layer.  The  dimensions  of 
this  layer  cannot  be  defined  precisely.  However,  a  generally  accepted 
practice  is  to  consider  the  boundary  layer  to  be  a  few  Debye  lengths 
thick.  The  Debye  length  is  the  characteristic  distance  of  a  plasma 
defined  by  the  relaxation  length  of  exponential  screening.  That  is,  it 
is  the  distance  that  is  required  for  the  electric  potential  to  drop  by  a 
factor  of  1/e.  The  significance  of  the  Debye  length  lies  in  the  fact 
that  the  electrostatic  potential  felt  by  an  electron  at  a  distance  much 
greater  than  a  few  Debye  lengths  from  a  charge  distribution  will  be  very 
small.  Because  the  Debye  length  is  the  characteristic  distance  of  a 
plasma,  it  is  natural  to  use  it  as  the  defining  parameter  for  the 
boundary  layer. 
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The  boundary  layer  sets  up  a  potential  screen  between  the  electrons 
leaving  the  surface  and  those  which  are  oeyond  the  boundary  layer.  This 
screen  reduces  the  energies  of  electrons.  It  also  determines  electron 
flow  on  the  surface  of  the  satellite  by  affecting  the  electromagnetic 
fields  at  the  surface.  Some  of  the  surface  current  may  get  into  the 
satellite's  interior.  This  is  the  important  physical  quantity  of  design 
interest.  Therefore,  a  detailed  knowledge  of  the  surface  currents, 
hence  boundary  layer,  is  an  important  step  to  the  understanding  of  the 
SGEMP. 

PRIOR  RELATED  WORK.  Serious  technical  effort  to  understand  the 
fundamental  physical  processes  of  SGEMP,  including  the  boundary  layer 
and  surface  currents,  began  in  the  mid  1960's.  The  Air  Force  Ueapons 
Laboratory  and  the  Defense  Nuclear  Agency  sponsored  large  SGEMP  programs 
to  develop  theoretical  knowledge  and  computational  techniques.  Conrad 
Longmlre,  Neal  Carron  and  others  analyzed  the  composition  of  the 
electron  boundary  layer  via  theoretical  and  semi -analytical  means  (Ref 
3,4,5,  and  6).  Others  concentrated  on  the  development  of  one-,  two-, 
and  three-dimensional  computer  codes  (Ref  7,8,9,  and  10).  An  excellent 
article  by  Higgins  et  al  (Ref  11)  reviews  the  progress  made  in  SGEMP 
analysis  up  to  1978. 

As  Higgins  points  out,  the  majority  of  numerical  approaches  used  to 
date  are  particle  simulation  codes.  This  technique,  also  called 
"particle  pushing",  relies  upon  a  statistical  model  of  macro-charges. 
That  is,  one  makes  the  assumption  that  a  given  amount  of  charge  can  be 
modelled  as  one  large  macro-particle.  This  allows  the  force  on  each 
particle  to  determine  its  movement  about  the  satellite.  Each  individual 
particle  represents  millions  of  electrons.  The  greater  the  number  of 
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macro-particles  used,  the  better  the  accuracy. 

Only  a  few  researchers  have  approached  transient  SGEMP  using  a 
continuum  model  for  the  electron  distribution  (Ref  12  and  13).  This 
approach  is  very  common  in  plasma  studies,  where  the  electrons  exhibit 
periodic  motion  (Ref  14  and  15).  In  this  method,  a  distribution 
function  is  defined  which  depends  upon  the  position,  x,  and  velocity,  v, 
of  the  electrons.  This  distribution  function  then  describes  the  density 
of  the  electrons  in  the  (x,v)  phase  space. 

The  major  difficulties  associated  with  the  particle-pushing  codes 
are  noise  fluctuations  produced  by  the  statistical  nature  of  the  algo¬ 
rithms,  and  the  computer  memory  requirements.  The  distribution  function 
approach  also  has  storage  problems.  Additionally,  the  dominant 
advection  terms  in  the  equations  are  difficult  to  treat  numerically. 
However,  as  Holland  points  out  (Ref  12),  for  cases  that  require  tracking 
of  many  particles  which  do  not  leave  the  numerical  mesh,  the  distribu¬ 
tion  function  approach  will  become  more  efficient  than  particle  pushing 
after  relatively  few  time  cycles.  Holland  roughly  estimates  that,  for  a 
particular  two-dimensional  problem  he  considered,  80,000  storage 
locations  would  be  required  by  a  particle-pushing  method  to  achieve 
results  of  comparable  accuracy  to  those  obtained  by  a  distribution 
tunction  approach  with  35,000  locations.  Thus,  the  distribution 
function  method  may  be  more  efficient  for  long  running  problems  with 
large  boundary  layer  effects.  But  the  advective  nature  of  the  equation 
remains  a  problem. 

THE  FINITE  ELEMENT  METHOD  (FEM).  The  FEM  is  a  numerical  technique 
which  has  been  used  extensively  by  engineers  to  perform  complex 
structural  analyses  in  the  aircraft  industry.  It  was  first  introduced 
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for  this  purpose  twenty  years  ago  (Ref  16).  Since  that  time,  research¬ 
ers  in  other  disciplines  have  recognized  that  for  some  problems  the  FEM 
has  certain  desirable  qualities.  As  a  result,  the  technical  literature 
has  become  filled  with  finite  element  applications  for  such  diverse 
fields  as  nuclear  reactor  technology,  electromagnetic  field  scattering 
and  oceanic  water  wave  behavior.  This  widespread  interest  in  the  FEM 
prompted  mathematicians  to  study  the  technique  in  great  detail.  Conse¬ 
quently,  there  are  now  many  texts  which  describe  the  FEM  from  both  the 
engineering  and  mathematical  points-of-view.  References  16,17,  and  18 
are  just  a  few  examples  of  these  texts. 

The  success  of  the  method  was  clouded  by  one  fact:  although  finite 
elements  worked  well  for  elliptic  and  parabolic  equations,  there  was 
some  doubt  as  to  the  usefulness  of  the  technique  on  hyperbolic  equations 
(Ref  17,  Chapter  7)  and  advective-dominated  transport  equations.  The 
equations  of  SGEMP  formulated  with  a  continuous  distribution  function 
for  the  electrons'  behavior  are  dominated  by  advection  terms.  Recently, 
several  authors  have  shown  that  the  FEM  can  be  effective  for  hyperbolic 
equations  and  for  advective-dominated  propagation  such  as  found  in  SGEMP 
(Ref  19,  Chapters  2  and  19). 

There  has  also  been  a  continuing  debate  over  the  accuracy  of  the 
FEM  versus  a  finite  difference  solution  to  the  same  equations.  However, 
Gresho  et  al  (Ref  19,  Chapter  19)  have  reported  improvements  over  finite 
difference  results  for  pure  advection  in  one  dimension  when  using  finite 
elements.  Also,  Demerdash  and  Nehl  (Ref  20)  claim  to  have  achieved  more 
accurate  results  for  static  and  sinusoidally  varying  nonlinear  electro¬ 
magnetic  field  problems  using  less  computer  time  and  memory  with  the 
FEM.  These  last  achievements  are  important  for  SGEMP  research  because 
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it,  too,  is  a  nonlinear  electromagnetic  problem. 

PURPOSE  OF  DISSERTATION.  The  major  goal  for  this  dissertation  is 
to  investigate  the  FEM  as  a  new  approach  for  the  one-dimensional  SGEMP 
problem.  The  weight  and  trial  functions  for  the  FEM  will  be  developed 
with  attention  given  to  the  strongly-advective  equations  of  SGEMP.  The 
elements  which  reduce  implementation  difficulties  will  be  used.  Param¬ 
eters  introduced  when  using  advective  methods  will  be  analyzed  to  de¬ 
termine  which  values  to  use.  The  nonlinear  SGEMP  equations  will  be 
integrated  using  methods  proven  to  work  for  other  finite  element  equa- 
t ions. 

The  impact  of  the  FEM  on  the  solutions  to  specific  SGEMP  prob¬ 
lems  will  be  examined.  The  sensitivity  of  the  finite  element  solutions 
to  the  space,  velocity  and  time  intervals  shall  be  analyzed.  Addition¬ 
ally,  results  using  this  different  technique  will  be  compared  with  the 
more  traditional  methods  of  analysis;  namely,  theoretical  solutions 
and  finite  difference  methods.  Thus,  another  set  of  results  will  be 
produced  to  describe  the  effects  of  the  nonlinear  boundary  layer. 

OVERVIEW.  The  finite  element  solution  to  the  one-dimensional  SGEMP 
equations  was  done  on  a  Control  Data  CYBER  176  series  computer  using  a 
FORTRAN  IV  program  called  FEMNEP  (_Finite  ^Element  Method  for  a 
^lonlinear  Electromagnetic  _Problem).  This  program  was  written  by 
the  author  as  a  research  tool,  not  specifically  designed  for 
large-scale,  production  running. 

Chapter  II  of  this  dissertation  presents  the  equations  and  physics 
of  SGEMP.  Prior  solution  techniques,  as  well  as  the  inherent  difficul¬ 
ties  built  into  the  equations  are  discussed.  Chapter  III  contains  the 
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basic  numerical  prescriptions  used  in  this  dissertation  to  solve  the 
equations  of  Chapter  II.  The  FEM  itself,  the  time  integration,  and  the 
special  numerical  treatments  are  included  in  this  chapter.  Chapter  IV 
presents  the  results  of  the  FEM  when  applied  to  the  one-dimensional 
SCEMP  boundary  layer.  Finally,  Chapter  V  summarizes  the  results  and 
lists  the  conclusions  of  this  dissertation. 

Many  of  the  derivations,  and  other  discussions  not  directly  related 
to  the  development  of  the  subject  matter  being  presented,  are  relegated 
to  the  appendices. 
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II.  System  Generated  Electromagnetic  Pulse 


THE  EQUATIONS.  For  the  purpose  of  this  study,  only  SGEMP  produced 
by  X-rays  will  be  considered.  Thus,  consider  an  object  illuminated  by  a 
pulse  of  X-rays  of  short  duration  (  ~1  psec)  in  a  vacuum.  Through 
nuclear  interaction  processes  with  the  object,  primarily  the  photoelec¬ 
tric  effect  for  Aluminum  (Ref  2  Chapter  25),  electrons  will  be  ejected 
off  its  surface.  The  behavior  of  these  electrons  in  their  own  self 
fields  is  determined  by  the  Vlasov  Equation,  a  form  of  the  Boltzmann 
Equation  (Ref  12  and  Ref  22,  p.  260): 


Since  the  electrons  of  interest  will  be  traveling  at  less  than  1% 
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of  the  speed  of  light,  the  acceleration,  a,  is  given  by 
Force, 


a(x,v,t)  =  -  £[E(x,t)  +  v  x  B(x, t ) ]  (2) 

ra 


for  rationalized  MKS  units.  In  this  equation. 


-1  9 

e  =  charge  of  electron  =  1.6  x  10  Coulomb 

-3  1 

m  =  mass  of  electron  =  9.11  x  10  Kgin 

£(x,t)  =  electric  field  in  Volts/m 

v  =  velocity  of  electron  in  i„,  sec 

5(x,t)  =  magnetic  induction  field  in  Weber/m2 


In  order  to  complete  the  solution  of  Eqn  (1),  Maxwell 
in  free  space  are  needed: 


-r 

v*£(x,t)  =  p(x,t }  (3a);  v4(x,t)  =  0  (3b) 


V  x  jTtf.t)  =  -  ||(x,t)  (3c) 


V  X  £(x,t)  =  |j  J(x,  t )  +  u  e  ^<x,t)  (3d) 

O  O  On 


p(x,t)  =  charge  density  in  Coul/m^ 
l(x,t)  =  current  density  in  Amp/m2 


the  Lorentz 


s  Equations 
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U  =  permeability  of  free  space 
0  =  4tt  x  10-7  Henry/m 

=  permittivity  of  free  space 

=  8.85  x  10-12  Farad/m 


The  solution  to  Eqns  (3a)  and  (3b)  is  facilitated  by  the  fact  that  if  £ 
and  2  fields  are  found  which  satisfy  these  equations  at  any  time,  then 
the  advancement  of  these  solutions  in  time  via  Eqns  (3c)  and  (3d)  will 
produce  fields  which  always  satisfy  them  (Ref  23  Sec.  1.2). 

One  other  auxiliary  equation  is  required  to  complete  the  physical 
description  of  the  SGEMP  process;  that  is,  the  connection  between  the 
distribution  function,  f(x,v,t),  and  the  driving  term  of  Maxwell's 
Equations,  3(x,t),  the  current  density.  This  is  really  a  matter  of  def¬ 
inition,  and  is: 

+<» 


3(x,t) 


-e 


/ 


vf (x,v,t) 


dv 


(4) 


Its  more  recognizable  form  is  3=-nev,  where  n  =  number  of  electrons  per 
cubic  meter.  An  actual  integration  can  be  performed  by  using  upper  and 
lower  limits  which  are  chosen  so  that  f(x,v,t)  is  negligible  outside 
these  limits.  In  principle,  then,  Eqns  (1),  (2),  (3c),  (3d),  and  (4) 
provide  the  physics  of  SGEMP.  All  that  is  needed  to  complete  the 
solution  are  initial  conditions,  boundary  conditions,  the  specification 
of  the  source  characteristics  and  a  technique  to  solve  the  equations. 

One  numerical  method  of  solution  would  proceed  the  following  way: 
assume  all  electromagnetic  fields  are  zero  to  start  with;  then,  with  a 
knowledge  of  S(x, v,0+At ) ,  the  distribution  function,  f(x,v,0)  is 
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advanced  to  f(x,v,0+At)  by  Eqn  (1).  From  this  information,  the  source 
term  for  Maxwell's  Equations,  J(x,0+At),  is  specified  through  Eqn  (4). 
Then,  Eqns  (3c)  and  (3d)  are  solved  for  B(x,0+At)  and  E(x,0+At)  in  some 
leap-frog  fashion.  Once  E  and  B  are  determined,  the  acceleration,  a,  is 
known  at  advanced  time,  t=At,  from  Eqn  (2).  The  procedure  now  can  be 
repeated  for  t=t+2At,  and  so  on. 

Of  course,  there  are  many  subtleties  in  this  prescription  which 
make  an  actual  numerical  solution  extremely  difficult.  The  numerical 
solution  to  Maxwell's  Equations,  alone,  is  a  challenging  task  (Ref  24). 

ONE-DIMENSIONAL  APPROXIMATION.  In  order  to  study  the  properties  of 
the  highly  space-charge-limited  region,  or  boundary  layer,  of  the  SGEMP 
problem,  it  is  usually  only  necessary  to  consider  the  motion  of  the 
electrons  in  one  spatial  dimension.  That  is,  a  one-dimensional  approx¬ 
imation  is  valid  for  an  infinite  flat  emitting  surface.  And,  the 
surface  may  be  considered  flat  if  the  boundary  layer  thickness  is  small 
compared  to  the  linear  dimensions  and  radii  of  curvature  for  the  region 


of 

interest. 

Now,  a  typical  boundary 

layer 

for  an  aluminum 

surface 

and 

typical 

incident  I  keV  blackbody  > 

C-ray 

spectrum  with 

a  0.001 

O 

cal /(cm  -nanosec)  flux  has  a  boundary  layer  thickness  of  about  1  mil- 


limeter. 

Therefore,  it 

is  reasonable 

to  make  a 

one-dimensional 

model  in 

order  to  study 

tin.  boundary 

layer 

since. 

satellites  have 

many  exterior  linear  dimensions  which  are 

much 

larger 

than  millimeter 

size. 

In  one  dimension,  several  simplifications  to  Eqns  (1)  thru  (4)  can 
be  made  since  all  dependent  variables  are  functions  of  one  spatial 
direction  and  one  velocity  direction  only.  Let  these  variables  be  z  and 
v2=v.  It  is  easy  to  show  from  both  physical  and  mathematical  principles 
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that  no  magnetic  fields  can  exist  in  these  equations  when  the  problem  is 

limited  to  one  spatial  dimension  (see  Appendix  A).  Letting  E  =E,  Eqns 

z 

(1)  and  (2)  can  be  combined  to  give: 


|-^(z,v,t)  +  v|£(z,v,t)  -  -E(z,t)|i-(z,v,t)  =  S(z,v,t) 
a  t  3  z  m  3  v 


(5) 


Likewise,  if  J  -J,  Eqns  (3d)  and  (4)  vield: 
z 


—( 

3t 


Equation  (5)  is  the  one-dimensional  Vlasov  Equation  and  Eqn  (6) 
says  that  the  Maxwell  displacement  current  is  proportional  to  the 
ordinary  current  when  no  magnetic  field  exists  (Ampere's  Law  minus  a 
magnetic  field  but  with  displacement  current).  The  Vlasov  Equation  can 
be  thought  of  as  a  statement  of  particle  conservation  in  a 
two-dimensional  phase  space.  The  simplified  Ampere's  Law  describes  the 
build-up  of  the  electric  field  as  a  function  of  the  electron  current. 

There  is  an  alternative  to  using  Ampere's  Law  in  the 
one-dimensional  approximation.  Gauss'  Law,  Eqn  (3a),  can  be  used  to 
determine  the  electric  field.  That  is,  a  complete  description  for  E  can 
be  obtained  from: 


||(z,t)  =  ip(z,t)  (7) 
o 

This  is  possible  because  E  responds  only  to  changes  in  the  charge 
properties  of  density  and  current,  not  magnetic  fields,  in  one 
dimension.  The  charge  density  p(z,t)  can  be  determined  from  f(z,v,t) 
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through  the  relationship. 


Therefr  *,  an  equivalent  approach  is  to  use  Eqn  (5)  with  Eqns  (7)  and 
(8)  to  solve  the  one  dimensional  SGEMP  problem.  For  higher 
two-dimensional  and  three-dimensional  situations,  Ampere's  Law,  Eqn 
(3d),  must  be  used  because  magnetic  fields  are  present. 

SOLUTION  TECHNIQUES  FOR  ONE-DIMENSIONAL  SGEMP.  There  are  two  basic 
approaches  that  can  be  used  to  attack  SGEMP,  numerical  simulation  and 
numerical  solution.  A  third  possible  technique,  exact  analytic 
solutions,  is  only  feasible  for  very  restricted,  idealized  cases.  This 
is  true  because  even  in  one  dimension,  the  equations  are  a  set  of 
coupled,  nonlinear,  integro-dif ferential  equations  which  have  no 
closed-form,  analytic  solution  for  an  arbitrary  source.  Nonetheless, 
there  are  some  interesting  cases  which  do  have  analytic  or  semi -analytic 
answers  (Ref  4). 

Numerical  simulation  is,  by  far,  the  most  common  approach  taken  to 
solve  the  SGEMP  problem.  This  technique  is  equivalent  to  the  solution 
of  Eqn  (5)  by  the  method  of  characteristics  (Ref  25,  Chapter  8).  The 
characteristic  lines  for  the  Vlasov  Equation,  by  itself,  are  those  which 
exist  in  (t,z,v)  space  such  that  their  directions  are  determined  by  dt  = 
k,  dz  =  kv,  and  dv  =  ka,  where  k  is  a  constant  and  a  =  acceleration. 
Thus,  if  numerical  solutions  are  sought  such  that  v  =  dz/dt  and  a  = 
dv/dt,  the  characteristics  of  the  Vlasov  Equation  are  being  followed. 
Mote  that : 
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a 


dv 

dt 


-  =  -  -E(z,t)  (9) 

m  m 


is  a  representation  of  Newton's  Second  Law.  The  idea  in  numerical 
simulation  is  that  the  particles  (electrons)  move  in  accordance  with  Eqn 
(9)  and  the  definition  of  the  particles  position  as  a  function  of  its 
velocity,  v  =  dz/dt.  Therefore,  "particles”  representing  a  large  amount 
of  charge  can  be  moved  about  with  Eqn  (9),  and  the  electric  field  found 
with  Gauss'  Law  or  Ampere's  Law.  This  is  called  "particle  pushing". 

There  are  several  documented  computer  codes  that  employ  particle 
pushing  (Refs  7,  8,  10  and  15).  The  major  drawback  of  this  method  is 
that  the  smoothness  of  the  solution  depends  upon  the  statistical  number 
^  of  "particles"  which  are  tracked  at  one  time.  This  is  why  it  is  truly 

numerical  simulation.  Electrons  are  simulated  by  a  given  amount  of 

charge  in  a  numerical  mesh.  They  are  injected  into  the  mesh,  moved 

I 

about  using  prescribed  physical  laws,  and  "killed”  when  they  leave  the 
system,  return  to  the  emitting  surface,  or  are  reduced  in  energy  below 
|  some  pre-determined  value.  The  primary  concern  in  actually  implementing 

such  a  procedure  is  the  particle  weighting  scheme;  that  is,  proper 
treatment  of  the  particle  density  within  each  cell.  The  actual  solution 
of  the  equations  of  motion  and  Gauss'  Law  (or  Ampere's  Law)  is  done 
|  through  conventional  finite  difference  techniques.  Elaborate  computer 

programs  have  been  written  which  solve  the  SGEMP  problem  in  one,  two, 
and  even  three  spatial  dimensions  using  particle  pushing  algorithms. 

■  An  alternative  to  particle  simulation  for  the  general  problem  is  to 

i 

!  solve  the  equations  numerically.  There  has  been  one  documented  finite 

difference  solution  to  the  plasma  equations  (Ref  12).  The  solution  was 
I  done  in  two  spatial  dimensions,  primarily  to  study  secondary  electron 

I 


I 
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effects,  not  highly  space-charge-limited  regions,  although  it  clearly 
could  be  used  for  such  a  purpose.  There  does  not  appear  to  be  a  well- 
documented  finite  difference  study  of  a  one-dimensional  SGEMP  model 
using  the  Vlasov  Equation  to  analyze  the  transient  behavior  of  the  boun¬ 
dary  layer. 

The  subject  of  this  dissertation  is  the  application  of  a  third 

numerical  technique,  the  Finite  Element  Method,  on  the  highly 

space-charge-limited  SGEMP  problem  in  one  dimension. 

MATHEMATICAL  CONSIDERATIONS.  Equation  (5)  with  Eqn  (6)  is  a  system 

of  nonlinear  integro-dif f erent ial  equations.  Certain  features  of  these 

equations  must  be  kept  in  mind  before  any  attempt  is  made  to  solve  them. 

The  Vlasov  Equation  describes  convective  transport  and  contains  the  ad- 
3  f 

vection  term,  v-r—  .  It  is  well  known  that  the  numerical  treatment  of 

a  Z 

such  a  term  requires  care.  For  example,  Richtmyer  and  Morton  show 
that  the  standard,  centered,  second  order  finite  difference  approxima¬ 
tion  of  this  term  is  not  satisfactory  (Ref  26,  Sec.  12.3). 

Also,  if  an  explicit  finite  difference  scheme  is  used  on  the  time 
variable,  t,  then  it  is  likely  that  some  type  of  restriction  must  apply 
to  the  relative  size  of  the  z-step  versus  t-step,  such  as  the  Courant- 
Friedrichs-Lewy  condition  (Ref  26  Sec.  10.2): 


-VAt  <;  , 
Az  1 


(10) 


This  condition  ensures  that  the  propagating  pulse  is  not  allowed  to 
traverse  more  than  one  z-cell  in  one  time  step  when  an  explicit  scheme 
is  used.  It  seems  reasonable  to  assume  that  a  similar  condition  should 


16 


exist  in  the  velocity  direction,  such  as: 


(11) 


where  a  is  the  acceleration.  However,  there  is  no  guarantee  that  any  of 
these  criteria  will  preserve  stability  for  the  finite  element  equations 
needed  in  SGEMP  calculations.  In  Chapter  IV,  data  will  be  presented 
which  demonstrate  the  usefulness  of  these  Courant  conditions  for  the 
present  problem.  There  are  some  cases  for  which  these  conditions  are 
unimportant.  This  will  also  be  discussed  in  further  detail  in  Chap¬ 
ter  IV. 

PHYSICAL  CONSIDERATIONS.  Besides  the  purely  mathematical  concerns 
mentioned  above,  there  are  physical  restrictions  which  need  to  be  con¬ 
sidered.  Since  the  distribution  function,  f(z,v,t),  represents  a 
particle  density,  physically  it  must  always  be  positive.  However,  when 
numerical  techniques  are  used  to  propagate  a  pulse  through  a  mesh,  os¬ 
cillations  in  front  of,  and  behind  the  pulse  often  occur  (Ref  27).  Con¬ 
sequently,  the  distribution  function  may  become  negative.  However,  this 
is  due  to  the  numerical  behavior  of  the  algorithms  rather  than  some 
actual  physical  effect.  It  has  been  reported  that  in  nonlinear  problems 
such  as  this  one,  the  negative  oscillations  can  very  quickly  become 
intolerable  (Ref  22,  p.  369).  Therefore,  steps  must  be  taken  to 
mitigate  this  potential  hazard. 

The  second  physical  requirement  on  -n  solution  to  the  SGEMP 
equations  is  due  to  the  concept  of  the  plasma  Debye  length.  In  ration¬ 
alized  MKS  units,  the  Debye  length  for  the  electrons  of  a  plasma  is 
defined  as  (Ref  28,  Chapter  8): 
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where  k.  =  Boltzmann  constant  =  1.38  x  10-23  Joule/°K 
T  »  temperature  of  electrons  in  °K 
n  =  number  density  of  electrons  in  m-3 


This  definition  is  based  upon  the  plasma  being  in  Maxwell-Boltzmann 
equilibrium,  and  represents  the  distance  within  which  the  electrons  are 
able  to  interact  individually.  The  number  of  electrons  within  a  Debye 
length  of  each  other  are  just  those  contained  within  a  sphere  of  radius 
=  X  .  Using  Eqn  (12),  this  turns  out  to  be  N  ,  where, 

D  A 


N 


4n  (Epk-T/e)372 


(13) 


The  charge  at  the  center  of  this  sphere  is  screened  from  those  which  lie 
outside  the  sphere. 

This  collective  behavior  of  a  plasma  has  important  implications  to 
the  numerical  solution  of  Vlasov's  Equation.  If  the  spatial  grid  is 
chosen  so  that  it  is  larger  than  the  Debye  length,  the  numerical  scheme 
will  not  be  able  to  handle  the  steep  gradients  in  particle  density  which 
would  occur  for  every  cell  width.  Therefore,  the  mesh  spacing  of  the 
spatial  variable,  Az ,  must  be  sufficiently  fine  to  resolve  the  plasma 


Debye  length 


BOUNDARY  AND  INITIAL  CONDITIONS.  In  order  to  complete  the  solution 
to  Eqns  (5)  and  (6),  boundary  conditions  and  initial  conditions  must  be 
presribed  for  the  system  of  equations.  The  initial  conditions  for  the 
problem  specify  that  the  following  variables  are  zero  at  t  =  0: 


f(z,v,0)  =  E(z,0)  =  S(z,v,0)  =  0,  for  all  v,z  >  0  ( 1 A ) 

The  boundary  conditions  on  the  Vlasov  Equation  can  be  treated  in 

the  same  manner  that  any  transport  phenomenon  is  handled.  In  the  (z,v) 

mesh  for  the  one-dimensional  problem  of  Fig.  1,  there  is  only  one 

physical  boundary  that  represents  the  division  between  the  infinite 

sheet  of  material  and  space,  the  boundary  at  z=0,  v  >  0.  However,  all 

four  boundaries  must  be  considered.  Basically,  free  surface  boundaries 

exist  for  this  equation  on  all  sides  of  the  mesh,  similar  to  the 

conditions  described  by  Bell  and  Glasstone  (Ref  29,  Sec.  l.ld).  The 

difference  is  that  for  the  Vlasov  Equation,  particles  are  restricted  to 

move  in  only  one  direction  for  any  given  part  of  the  mesh.  Therefore, 

particles  which  are  allowed  to  enter  (exit)  the  mesh  at  z=0  (  z=z  ) 

max 

cannot  exit  (enter)  from  that  same  part  of  the  boundary.  This  makes  the 

free  surface  boundary  conditions  exactly  correct  at  these  locations,  not 

the  idealization  referred  to  by  Bell  and  Glasstone.  Using  the  free 

surface  boundary  conditions  at  z=±v  is  an  idealization,  however,  for 

max 

the  assumption  must  be  made  that  no  electrons  can  return  after  leaving 
the  mesh  from  either  of  these  interfaces.  Specifically,  the  free 
surface  boundary  conditions  require  that  electrons  which  leave  the  mesh 
cannot  return.  Therefore,  the  following  conditions  must  be  met  along 
the  boundaries  of  the  mesh,  as  depicted  in  Fig.  1: 
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f(z,+v  ,t) 
max 


f(z  ,v,t)  = 
max 


f (0 ,v, t )  = 


f(z,-v  ,t) 
max 


0,  for  all  z,t  (15a) 

[  0,  v<0 

free,  vX) 
source,  vX) 
free.  v*0 
free,  for  all  z,t  (15d) 


(15b) 


(15c) 


Equation  (15a)  is  true  because  it  is  assumed  that  no  electrons  can 

accelerate  out  the  top  of  the  mesh.  Therefore,  the  only  flow  possible 

is  into  the  mesh,  which  is  set  to  zero.  The  distribution  function  is 

allowed  to  be  whatever  it  wants  to  be  at  z  ,  v  >0,  and,  z=0,  v  *  0; 

max 

that  is  a  free  boundary. 

Special  note  must  be  made  about  the  conditions  that  exist  at  z=0, 
vX>,  Eqn  (15c).  Because  the  source  of  electrons  is  assumed  to  be  only 
at  the  surface,  >_he  boundary  condition  at  this  interface  is  the  true  re¬ 
quirement  to  get  electrons  into  the  mesh.  That  is,  there  is  no  volume 
source,  S(z,v,t),  which  is  spread  throughout  the  region  of  interest. 
Rather,  the  source  only  exists  at  the  boundary,  and,  therefore,  actual 
insertion  of  electrons  into  the  mesh  will  be  done  through  a  boundary 
s  pecif ication. 

The  boundary  condition  for  the  electric  field  is  taken  care  of 
automatically  by  having  the  correct  distribution  function.  There  are 
only  two  boundaries  for  E(z,t):  z=0  and  z=  <*>.  At  the  interface  of  a 
near  perfect  conductor  (like  Aluminum)  the  boundary  conditions  for 
Maxwell's  Equations  require  that  only  the  normal  component  of  E  exist 
near  the  surface.  Inside  the  conductor,  E  is  zero  (Ref  30,  Sec.  8.1). 
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Therefore,  the  E  field  on  the  surface  is  determined  by  having  the 

correct  charge  density  on  the  surface.  Since  it  is  reasonable  to  assume 

a  perfect  conductor  for  the  SGEMP  problem,  the  distribution  function, 

f(0,v,t),  a  surface  charge  density,  will  determine  the  correct  electric 

field  on  the  surface.  At  the  other  boundary,  z=z  ,  the  electric 

max 

field  is  allowed  to  take  on  whatever  value  f(z  ,v,t)  forces  it  to  be. 

max 

SOURCE  DESCRIPTION.  The  final  remaining  topic  of  discussion  to 
conclude  this  chapter  on  SCEMP  is  the  way  that  the  source  of  electrons 
can  be  described.  The  determination  of  the  electron  emission  properties 
is  a  very  complicated  process  in  itself.  There  are  elaborate  computer 
programs  which  perform  these  calculations  (Ref  31).  They  take  into 
account  such  factors  as:  the  thickness,  shape  and  composition  of  the 
target;  the  angle  of  incidence,  energy  spectrum,  time  history  and  total 
fluence  of  the  photon  radiation;  and  the  individual  nuclear  processes 
which  release  electrons  (Compton  scattering,  photoelectric  effect  and 
pair  production).  This  is  a  difficult  problem  to  solve  in  itself.  For¬ 
tunately,  however,  it  is  possible  to  separate  the  electron  emission 
problem  from  the  SGEMP  problem.  That  is,  the  electrons  which  are 
emitted  into  space  from  initial  arriving  radiation  do  not  interact  ap¬ 
preciably  with  the  late  arriving  radiation.  The  probability  of  a  photon 
interacting  with  an  electron  in  free  space  via  Compton  or  Thomson 
scattering  is  extremely  small  due  to  the  very  low  density  of  electrons. 
For  example,  the  electron  number  density  for  a  typical  SGEMP  plasma  is 
n  ~10‘  cm  .  For  Thomson  scattering,  the  cross  section  for  interac¬ 
tion  is  6.0  x  10~'-s  cm^/e  ,  which  is  larger  than  for  Compton  scattering 
(Ref  2,  Chapter  23).  Therefore,  the  linear  attenuation  coefficient  is 

—  |  o 

about  1.0  x  10  cm  at  most.  This  is  to  be  contrasted  with  photon 
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scattering  probabilities  in  a  solid  material  with  densities  of  I024  cm-3 , 
and  linear  attenuation  coefficients  of  approximately  0.1  cm-1  .  As 
a  result,  all  treatments  of  SGEMP  in  existence  rely  on  pre-determined 
emission  properties. 

As  was  just  noted  in  the  previous  section,  the  SGEMP  problem  has  a 
unique  condition  on  the  source  of  electrons.  Because  there  is  no  actual 
volume  source,  but  only  a  source  defined  at  an  interface,  the  emission 
of  electrons  is  really  a  boundary  condition.  Thus,  from  now  on,  Eqn  (5) 
will  be  considered  to  have  no  driving  term  on  the  right-hand-side. 
Rather,  the  boundary  conditions  on  Eqn  (5)  will  act  as  the  electron 
driver. 

An  assumption  that  is  very  frequently  made  about  the  electron 
source  is  that  the  energy,  hence  velocity,  dependence  and  time 
dependence  are  separable.  That  is,  one  assumes, 

f(0,v,t)  =  T(t)  fv(v),  v>0  (16) 

This  separation  is  for  convenience  only,  and  is  approximately  true. 
Should  some  variation  of  the  velocity  spectrum  as  a  function  of  time  be 
supplied,  there  are  no  unique  difficulties  which  would  arise  in  the 
solution. 

VELOCITY  DEPENDENCE  OF  ELECTRON  SOURCE.  Numerical  calculations  of 
the  photomission  properties  for  X-ray  radiation  show  that  the  electron 
energy  spectrum  is  represented  well  by  an  exponential  function  when  the 
X-rays  are  treated  as  blackbody  radiation.  The  results  also  show  that  a 
a  cos9  angular  distribution  is  a  good  approximation  (Ref  5).  If  the  X- 
rays  are  normally  incident  on  a  flat  surface,  then  the  electron  velocity 
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spectrum  is  reduced  to  the  following  expression  In  one  dimension  (see 
Appendix  B): 


f  (v) 
v 


[sec2 /m4  ] 


(17) 


where  V1  =»  material  yield  (electrons/Joule ) 

w^  =  exponentiation  energy  (Joules), 

which  is  property  of  material  and 
photon  spectrum 


4>0  =  X-ray  fluence  (Joule/m2) 

E1 (x)  =  exponential  integral  (Ref  32,  Chapter  5) 


00 


It  should  be  noted  that  variations  other  than  the  above  type  of  exponen¬ 
tial  dependence  of  the  velocity  spectrum  are  possible,  depending  on  the 
exact  nature  of  the  source. 

TIME  DEPENDENCE  OF  ELECTRON  SOURCE.  The  choice  of  the  time 
dependence  of  the  source  function,  T(t),  is  a  crucial  one  since  it  is 
this  parameter  which  determines  whether  or  not  scaling  of  the  SCEMP 
problem  is  possible.  Carron  and  Longmire  (Ref  6)  have  shown  that  if  the 
X-ray  pulse  rises  as  a  integer  power  of  time,  the  one-dimensional, 


normal  incidence  SGEMP  boundary  layer 

scales. 

This  means 

that 

the 

problem  can  be  solved  once,  with 

all 

variables 

scaled  to 

this 

one 

solution.  If  the  scaled  variables 

are 

represented 

by  t ' ,  z  ' , 

v. 

f, 

a',  and  S',  then  these  similarity  variables  are  defined  by:  t=Tt', 
z=Lz',  v»<v>v',  f-Nf'/<v>,  a«<v>a'/T,  and  S-NS'.  The  variables,  t,  z. 


24 


v,  f,  a  and  S  are  defined  in  Eqns  (2)  and  (5).  The  parameters,  T,  L, 
<v>  and  "  are  the  plasma  period,  characteristic  length,  average  velocity 
and  a  reference  number  density  of  the  electrons,  respectively.  If  the 
time  dependence  is  slowly  varying  enough,  then  the  equations  may  be 
reduced  to  a  quasi-static  case  fof  which  closed-form  or  quadrature 
solutions  exist  (Ref  A). 

Although  these  cases  are  of  general  interest  for  their  simplicity 
and  for  certain  laboratory  sources,  SGEMP  caused  by  a  nuclear  detonation 
and  some  pulsed  electron  beam  sources  are  driven  by  an  approximately 
exponentially  rising  and  falling  pulse.  One  analytic  form  of  this 
dependence  frequently  used  is: 


T(t)  =  c3 


_ e 

i  +  e 


cit 

=2  > 


(18) 


where 


c3 

C1 

c2 


normalization  constant  (sec-1) 
approximate  rise  rate  (sec-1  ) 
approximate  fall  rate  (sec-1  ) 
approximate  peak  time  (sec) 


The  following  condition  is  used  as  a  normalization  requirement: 


+<x> 

j T(t)  dt  =  1 

— ao 
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When  this  integral  is  carried  out,  c3  is  determined  to  be. 


a,  sin(7?c,  / c.^  ) 


If  the  source  is  represented  by  Eqns  (16),  (17)  and  (18),  only  a 
general  numerical  solution  is  possible  using  either  direct  techniques  or 
particle  simulation. 

SUMMARY  OF  NUMERICAL  PROBLEM.  The  next  chapter  begins  the  discus¬ 
sion  of  the  finite  element  techniques  for  the  SGEMP  boundary  layer  in 
one-dimension.  For  ease  of  reference,  the  mathematical  problem  is  re¬ 
stated  here.  The  distribution  function  for  the  electrons,  f(z,v,t),  sa¬ 
tisfies  the  Vlasov  Equation: 


i^(z,v,t)  +  vi!(z,v,t)  - -E(z,t)|-^(z,v,t)  =  0  (5) 

3t  3z  m  3v 


The  electric  field,  E(z,t),  is  related  to  f(z,v,t)  via: 


vf(z,v,t)  dv  (6) 


all  v 


The  initial  conditions  for  this  system  of  equations  are  that: 


f(z,v,0)  =  E(z,0)  *  0,  for  all  v,z  >  0  (14) 
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The  boundary  conditions  needed  are: 


f(z,+v  ,t)  =  0,  for  all  z,t  (15a) 

max  ’  ' 


f(z  ,v,t)  = 
max  ’ 


0,  v  <  0 
free,  v  >0 


(15b) 


f (0,v,t) 


I  T(t)f  (v), 

(  free,  v  <  ( 


(15c) 


f(z,-v  ,t)  =  free,  for  all  z,t  ( 1 5d ) 

max 


The  behavior  of  the  distribution  function  at  the  boundary  can  be  approxi¬ 


mated  by: 


(17> 


Finally,  the  time  history  of  the  electron  source,  T(t),  is  assumed  to  be 
of  the  form: 


T(t)  =  c  — °8) 
1  +  e  2  u 
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III.  The  Numerical  Techniques 


This  chapter  will  present  all  the  numerical  approaches  used  for  the 
solution  to  the  SGEMP  one-dimerssional  boundary  layer  problem  discussed 
in  Chapter  II.  Primarily,  this  consists  of  the  application  of  the 
finite  element  method  (FEM)  to  Eqns  (5)  and  (6).  But  there  are  several 
other  considerations.  The  time  integration  is  a  vital  part  to  the 
solution  of  this  problem.  Also,  the  boundary  conditions  must  be  handled 
with  care.  Finally,  the  resulting  set  of  algebraic  equations  requires 


an  efficient  solution 

technique , 

which  is  closely 

coupled 

to  the 

computer 

requirements 

(storage  and 

processing  time) 

for  the 

entire 

method. 

First  of  all, 

a  discussion 

of  the  FEM,  itself, 

is  in 

order. 

THE 

FINITE  ELEMENT  METHOD.  The 

FEM  is  a  numerical 

technique 

which 

provides 

a  systematic 

method  for 

solving  physical 

problems 

in  an 

arbitrary 

global  mesh. 

The  mesh 

can  be  nonuniform 

and  irregular, 

composed 

of  "elements" 

of  any  shape, 

but  usually  they 

are  chosen  as 

simple  polygons  (rectangles,  triangles,  tetrahedrons,  etc.).  Another 
important  feature  of  the  technique  is  that  it  allows  controlled  applica¬ 
tion  of  virtually  any  reasonable  order  polynomial  approximation  in  the 
mesh.  Mathematically,  the  FEM  is  an  extension  of  the  Ray- 
leigh-Ritz-Galerkin  technique  (Ref  17,  Chapter  1).  The  global  region  is 
subdivided  int-  smaller  regions  called  "elements".  The  unknown  function, 
g(x),  is  expanded  into  a  set  of  piecewise  polynomials ,  (x ) .  That  is, 


g(x) 


(21) 
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The  FEM  prescibes  the  way  that  the  trial  functions  are  to  be  chosen. 

For  Lagrangian  elements,  they  are  built  so  that  for  every  node,  k, 

')>  .(x,  )  =  5  .  This  gives  the  expansion  coefficients,  g  ,  physical 
J  k  jk  j 

significance.  In  this  dissertation,  the  expansion  coefficients  are  the 
distribution  function  values,  and  the  the  electric  field  values,  at  the 
nodes. 

METHOD  OF  WEIGHTED  RESIDUALS.  There  are,  in  general,  two  different 
ways  that  the  solution  via  the  FEM  for  any  problem  may  be  set  up  (Ref 
18,  Sec.  3.4).  The  first  approach  is  to  solve  a  variational  problem. 
That  is,  the  physical  situation  is  stated  as  an  integral  relationship 

for  an  entire  set  of  functions  -  a  functional.  Then,  the  correct 

solution  is  that  function  which  minimizes  the  functional.  This  method 
is  very  common  when  there  exists  some  physical  variable  which  must  be  a 
minimum  as  some  parameter  is  varied  throughout  its  range,  such  as  the 
Lagrangian  (kinetic  energy  -  potential  energy)  of  a  conservative  system. 
The  FEM  is  an  integral  (global)  solution  technique,  as  opposed  to  a 
discrete  (local)  method  such  as  finite  difference  methods.  The  latter 
approach  numerically  solves  the  differential  equations  for  the  physical 
s ituat ion. 

Of  course,  both  the  differential  and  integral  equations  represent 
the  same  physical  reality.  However,  it  is  not  always  possible  to 
determine  the  correct  functional  for  a  particular  problem.  When  this 
happens,  a  different  approach  must  be  taken.  The  FEM  can  still  be  used 
by  developing  global,  integral  expressions  from  the  differential 
equations.  The  Method  of  Weighted  Residuals  (MWR)  is  such  a  process. 
The  residual  is  the  difference  between  the  true  solution  and  the  approx¬ 
imate  one.  The  differential  equation  is  multiplied  by  a  set  of 
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functions,  called  weights,  and  the  new  equation  is  integrated  over  the 
entire  solution  domain.  This  converts  a  local,  differential  equation 
into  a  global,  integral  relationship.  In  this  work,  the  finite  element 
equations  are  derived  using  the  MWR,  as  discussed  in  Zienkiewlcz  (Ref 
18,  Sec.  3.4). 

GLOBAL  FINITE  ELEMENT  EQUATIONS.  Let  the  distribution  function, 
f(z,v,t),  be  approximated  by  a  set  of  N  trial  functions,  N.(z,v), 
i=l , . . . ,N.  That  is , 


f(z,v,t)  =y^N.(z,v)f  .(t)  (22) 

where  f.(t)  is  a  nodal  value  of  f;  f  ,(t)=f (z  .,v  ,,t).  Likewise,  let 

3  3  3  3 

N 

z 

E(z,t)  Vz)Ek(t)  (23) 
k=l 

where  M, (z)  are  a  set  of  N  trial  functions  for  the  electric  field,  and 
k  z 

E  (z)=E(z  ,t)  are  the  E-field  nodal  values.  These  approximations  for 

K  K 

E(z,t)  and  f(z,v,t)  are  used  in  Eqn  (5).  When  the  result  is  multiplied 
by  W  (z,v),  i=l,...,N,  the  MWR  weights,  and  summed  and  integrated  over 
the  entire  region,  the  following  equation  is  obtained: 


ttff 

i-1  j=l  z'v 


[W.(z,v)K.(z,v)~j(t)  +  vWi(z,v)||j(z,v)f  .(t) 


M  (z)E  (t)W.(z,v)^ij(z,v)  f  .  ( t )  ]  dvdz  =  0  (24) 

K  K  1  iJ  V  J 
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Recall  that  the  source  function,  S(z,v,t),  has  been  set  to  zero  because 
the  electrons  are  emitted  via  the  imposition  of  an  appropriate  boundary 
condition.  This  equation  can  be  written  in  a  more  condensed  matrix 
form: 


[A]NxN|d£(t|  +  [a_jNxN  |f(t)j'Xl-  [Av]jXNEJ(t>  |f(t)}* 


=  0  (25) 


with  the  help  of  the  following  definitions: 


(A)..  =  /  /w. (z,v)N.(z,v)  dvdz  (26a) 

iJ  J  J  1  J 

z  v 

(Az)ij  =  J jv (z,v)  dvdz  (26b) 


[(A  )  ] 
v  J 


.  .  =  - /> 1  (z)[  /w  (z,v>~j(z,v)  dv] 
lj  m  I  J  /  l  ov 


dz  (26c) 


Thus,  [A],  [A  ],  and  the  [A  ]  are  NxN  matrices,  and  { f ( t ) )  is  an  Nxl 
z  v 

matrix  (a  column  vector). 

Once  the  weight  and  trial  functions  are  chosen,  and  the  integrals 
in  Eqns  (26)  are  performed,  Eqn  (25)  can  be  solved  as  a  set  of  N  first 
order,  ordinary  differential  equations  in  time.  Applying  the  same 
approach  to  Ampere's  Law  in  one  dimension,  Eqn  (6),  one  finds: 


N  xN  | 
,  z  z 


N  xN  <  i  N  xl 
z  2  f(t)  2  =  0  (27) 
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where. 


(D) 


ij 


P  (z)M  (z)  dz 

i  j 


(28a) 


<V 


.  .  =  -  P. (z)[  |  vQ.(z,v)  dv] 

1J  ZoJ  1  J  J 


dz 


(28b) 


and  where  the  dot  above  the  E  means  differentiation  with  respect  to 

time.  The  P.(z),  i=l,...N  ,  are  the  weight  functions;  Q  (z,v), 

1  z  J 

,  are  the  trial  functions  for  f(z,v,t);  and  the  M.(z)  are  the 
z  J 

trial  functions  for  E(z,t). 

Note  that  in  this  general  formulation  of  the  finite  element 
equations,  the  nodal  points  are  not  necessarily  the  same  in  both 
equations.  That  is,  they  could  be  solved  in  two  separate  meshes. 

ELEMENT  SHAPE  AND  MESH  DESCRIPTION.  Before  proceeding  further,  it 
is  necessary  to  discuss  the  choice  of  element  and  mesh  configuration, 
and  the  weight  and  trial  functions  which  are  compatible  with  these 
choices.  The  simplest  element  shape  that  can  be  used  for  the  Vlasov 
Equation  is  a  rectangular  element  in  (z,v)  space.  There  are  many 
advantages  in  using  rectangles  in  this  problem,  and  some  disadvantages. 
In  order  to  demonstrate  the  methods  useful  in  the  application  of  the  FEM 
to  one-dimensional  SGEMP,  I  chose  rectangles.  For  a  detailed  discussion 
of  the  advantages/disadvantages  of  different  element  shapes  and  less 
retricted  meshes,  see  Appendix  C. 

A  rectangular  mesh  is  regular.  That  is,  for  linear  Lagrangian 
elements,  the  node  pattern  consists  of  the  set  of  (z^  ,  v  )  with 
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1=1,... ,N  ;  J=1,...,N  .  For  every  z  there  are  N  v  's  and  for  every 

v  z  I  v  J 

v  ,  there  are  N  z  's.  A  typical  rectangular  mesh  is  depicted  in  Fig. 
J  z  I 

2.  Of  course,  the  mesh  can,  in  general,  be  non-uniform  in  either  the  z 
or  v  direction  independently. 

THE  TRIAL  FUNCTIONS.  The  description  of  the  trial  functions 
(sometimes  called  approximating  functions  or  shape  functions)  can  most 
easily  be  described  in  the  local  coordinate  system  of  the  element.  For 
rectangles,  the  simplest  trial  functions  which  can  be  used  are  linear 
polynomials  with  nodes  at  each  of  the  four  corners,  as  shown  in  Fig.  3. 
These  trial  functions  are  members  of  the  Lagrange  family  of  polynomials. 
They  are  also  belong  to  a  family  of  polynomials  which  are  especially 
suited  for  the  FEM  because  of  their  simplicity.  Zienkiewicz  refers  to 
the  entire  sec  as  the  Serendipity  family  (Ref  18,  Sec.  7.5).  They 
satisfy  the  essential  requirement  of  the  finite  element  approximation: 
continuity  of  the  unknown  function  across  element  interfaces.  These 
linear  functions  are  defined  in  the  local  coordinate  system  by: 

N^C.n)  =  1  +  n.n)  (29) 

i  =  1,...,4 

where  the  transformation  from  the  global  to  the  local  system  is  made  by: 

^  =  — ( z  -  z  )  (30a);  p  =  -kv  -  v  )  (30b) 

a  c  be 

The  constants,  a,  b,  z  and  v  are  defined  in  Fig.  3.  Note  that  each 

c  c 
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of  these  functions  has  the  property  of  being  unity  at  one  node,  and  zero 
at  the  other  three  nodes,  with  a  linear  variation  in  between. 

Of  course,  the  variation  in  the  £  and  n  directions  are  independent 
of  one  another,  and  as  such,  can  be  treated  separately.  Thus,  let 


VO  =  1(1  -  O  (31a) 

Lo(£)  =  ^(1  +  O  (31b) 

with  like  definitions  in  the  n  -direction.  These  are  the  basic  linear 
trial  functions  in  one  dimension,  sometimes  called  tent  functions 
because  they  construct  a  tent-like  shape  when  they  are  built  on  a  series 
of  nodes  (see  Fig.  4a).  The  rectangular  functions  of  Eqn  (29)  are 
built  from  these.  Using  the  nodal  numbering  system  of  Fig.  3: 


Nj(£>n)  =  lq(£)L,j(r\) 

(32a) 

N2(c,n)  =  VOVn) 

(32b) 

Nj  ( £  >  n )  =  ^(^^(n) 

(32c) 

Mc.h)  =  L1(c)L2(n) 

(32d) 

Just  as  the  combination  of  tent  functions  over  one-dimensional  elements 
produce  the  global,  piecewise  trial  functions  shown  in  Fig.  4a,  the 
global  functions  required  for  the  two-dimensional  problem,  N  (z,v), 

i»l,...N  are  built  from  the  local  trial  functions,  N^(£,n  )>  i=l . 4. 

Figure  4b  shows  the  tent  functions  for  one  element  in  two  dimensions, 
and  Fig.  4c  depicts  the  tent  functions  at  one  node  over  a  series  of 


rectangles. 

WEIGHT  FUNCTIONS  AND  UPWINDING.  The  residual  weights,  Wi(z,v), 
i=I,...N,  are  also  piecewise  polynomials  and  so  can  be  constructed  most 
easily  in  the  local  coordinate  representation.  One  of  the  more  common 
choices  for  the  weight  functions  is,  W^(z,v)  =  N^(z,v),  for  all  i.  This 
is  sometimes  referred  to  as  the  Galerkin  Method.  However,  an  upwinding 
procedure  will  be  used  (Ref  33,  and  Ref  19,  Chapter  1)  in  anticipation 
of  special  requirements  for  the  advection  term  of  the  Vlasov  Equation. 
In  this  approach,  Galerkin  weighting  appears  as  a  special  case. 
Referring  to  the  local  system,  let: 


wl<F,»n)  =  IM<0 
w2U,n)  =  [L2(0 
w3(5,n)  -  1^(0 

w„U,n)  =  (MO 


+  a1F(OJ[L1(n)  + 

-  ctiFUmLiCn)  + 

-  a3F(C)][L2(n)  - 
+  a3F(£)3 [L2(n)  - 


OqF(n)]  (33a) 

«2  F ( n ) ]  (33b) 

a2  F  ( n )  ]  (33c) 

cft*F(n)]  (33d) 


The  a^  are  the  upwinding  parameters.  Note  that  Galerkin  weighting 

occurs  when  a  =0,  for  all  i*l,...,4.  The  side  for  each  a  is  shown  in 
i  i 

Fig.  5.  The  function  F(*)  is  chosen  so  that  F(*)=0  at  the  nodes.  This 
preserves  continuity  of  the  weight  functions  at  element  interfaces.  If 
F(.)  is  a  positive  function,  the  sign  of  a.  is  picked  to  match  it  with 
the  direction  of  flow  of  electrons.  Thus,  for  a3  and  a3 ,  the  sign  will 
be  positive  when  v>0,  and  negative  when  v<0.  For  a2  and  ,  the  sign  is 
always  negative  because  the  electric  field  decelerates  the  electrons 
during  the  times  of  interest.  Finally,  the  integral  of  F(  •  )  over  the 
element  should  be  unity.  A  function  which  fits  all  these  requirements 
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is : 


F(x)  =  -  |(1  -  x ) ( 1  +  x)  =  -  3L1(x)L2(x)  (34) 

Equations  (33)  have  been  constructed  from  the  following  one  dimen¬ 
sional  expressions  for  each  side  of  the  rectangle: 

W  (x)  =  L  (x)  +  (-l)1+1ctF(x)  (35) 

i  =  1,2 

The  shape  of  one  of  these  weight  functions,  Wj(x),  is  .depicted  in  Fig. 
6  for  several  a's.  It  is  the  biasing  of  the  weights  towards  one  node 
with  repect  to  the  other  node  that  gives  upwinding  an  advantage  for 
propagation  problems. 

ELEMENT  MATRICES.  The  actual  evaluation  of  Eqns  (26)  and  (28)  is 
normally  done  in  the  local  coordinate  system  of  the  element.  Then,  the 
transformation  is  made  to  the  global  system.  This  transformation  is 
known  as  the  "assembly"  process  of  the  FEM  and  involves  a  Boolean 
mapping  of  the  local  numbering  system  for  the  nodes  to  the  global 
numbering  system  (Ref  34,  Chapter  6).  In  symbolic  form,  for  example, 

N 

(A]NxI'T  =^[a(e)]4x4  (36) 

e  =  l 

which  says  that  the  global  NxN  [A]  matrix  is  constructed  from  a  4x4 
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Fig.  6  Weight  Functions  in  ID 

For  Various  Values  of  Alpha 


3 


element  matrix  (for  linear  rectangular  elements)  by  "summing"  over  the 
number  of  elements,  N  ,  in  the  system.  The  script  -A"  is  used  to 
denote  that  this  sum  is  a  special  assembly  procedure  defined  by  a 
Boolean  transformation.  Appendix  D  has  the  details  on  the  exact  meaning 
of  Eqn  (36). 

The  expressions  for  [a^e^]  and  [a^]  are  straightfoward : 


+  1  +1 

(a(e))..=ab j  /w.(C,n)NX5,n)  dCdn  (37a) 

-1  -1 
+1  +1 

(a(e))..  =  b  f  /w.U.nXbn  +  V  >B  u,n)  d^dn  (37b) 

z  1J  J  J  1  c  3? 

-1  -1 


In  these  equations,  the  Jacobian  of  the  transformation  from  global  to 

local  coordinates,  a«b,  has  been  included. 

The  expression  for  matrix  [A  ]  ,  Eqn  (26c),  is  a  little  bit  more 

v  J 

difficult  to  write  in  local  coordinates  because  of  the  explicit 

dependence  on  each  M  (z).  However,  M  (z)  is  a  tent  function  over  the 

J  J 

Jth  node  for  the  linear  approximation  of  regularly-spaced  rectangles. 

Therefore,  M^(z)  in  Eqn  (29c)  is  either  Lj(£)  or  L 2  (  K  )  in  the  local 

system.  The  assembly  of  [A  ]  does  not  occur  through  Eqn  (36),  but 

v  J 

rather  via: 


[A  ]"/ 
V  J 


k'i 

=A< 


j(e )4xl*  f  a(e)<*4} 


That  is,  [A  j  is  built  from  two  element  matrices,  rather  than  one,  and 
v  J 

the  way  that  the  matrices  are  combined  depends  not  only  on  the  Boolean 


transformation,  but  also  on  the  Jth  position  of  the  global  column  of 
nodes.  See  Appendix  D  for  details. 

Explicitly,  then, 


+1  +1 


f a  ) 


i:j  =  */  /yc)w.(E,n)|ij(?,n)  d5dn 
_1  e  =  1,2 


Now,  using  Eqns  (32)  and  (33),  the  elemental  matrices  can  be 
evaluated  in  closed  form  to  obtain: 


(  a(e )  ■)  _  ab 

v  Jij  4 


(l+4^-^-+n-ct-^1+Tn-n-+^<ct-  i  ^  (40a) 

3  i  j  ii  3  i  j  i  l-l 


.a(e))  .  .  =*  2u  (1  +  5  a  )[(3  +  n  n .  +  3n.a  )v  +  (n  +  n.  +2n,  n.a  )b] 
z  ij  12  J  ik  i  j  ipe  l  J  5  i  j  p 

(40c) 


•gnj<(2 wy  <5+75i<9+^hi  <4°') 


,  ,  , .  ( 1  +  n  a  )  i 

^u'ro-j  <2  15  +  3h«;V 


(40d) 


where  i,j*l,...,4;  and  all  subscripts  are  cyclic  (that  is,  j— 1=4  when 
j=l,  and  j+2=i  when  j=3,  etc.).  The  following  replacement  rules  are  in 
effect  for  k  and  p: 

k  =  1  when  i  =  1 ,  2 
k  =  3  when  i  =  3,  4 
p  =  4  when  i  =  1 ,  4 
p  =  2  when  i  =  2 ,  3 
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Example  derivations  for  these  equations  are  shown  in  Appendix  E. 

CURRENT  DENSITY  CALCULATION.  The  only  remaining  element  matrices 

to  evaluate  are  [d^e^]  and  [d^e^]  which  are  associated  with  Ampere's 

Law.  There  is  a  process  that  can  be  used  for  a  rectangular  mesh  which 

eliminates  the  need  to  carry  out  the  integrals  of  Eqns  (28).  Everywhere 

that  a  nodal  value  of  E(z  ,t)  is  needed,  nodal  values  of  f(z  ,v  ,t) 

J  J  J 

exist.  Therefore,  a  simple  numerical  integration  in  the  v-direction  can 
be  used  to  evaluate  Eqn  (6)  which  will  determine  E(z^,t)  for  all  J.  In 
other  types  of  meshes,  this  would  not  be  possible  (see  Appendix  C). 

Since  f(z,v,t)  is  assumed  linear  in  v,  the  integration  required  in 
Eqn  (6)  is  quadratic  in  v.  Therefore,  Simpson's  rule  (Ref  35,  Chapter 
1A)  will  evaluate  the  integral  exactly  over  each  element  in  this  approx¬ 
imation.  Thus: 


E 

J 


(t) 


N  -1 
v 


v  )[2v  f(z  ,v  ,t)  +v  f(z  ,v 
k  k  J  k  k  J  k+1 


t) 


+  v  f(z  v  ,t)  +  2  v  f(z  ,v  , t ) ]  (Al) 

k+1  j’  k  k+1  J  k+1 

N  =  number  of  nodes  in  v-direction. 
v 

TIME  INTEGRATION.  Equations  (25)  and  (Al)  are  first  order, 
coupled,  ordinary  differential  equations  in  time.  The  time  integration 
is  performed  by  using  a  two-step  procedure  based  on  the  Lax-Wendroff 
method  (Ref  19,  Chapter  2).  The  algorithm  is  defined  by  the  following 
equations : 

f(t  +  41)  =  f(t)  +i!f(t)  (A2a) 

2  2 

f(t  +  At)  =  f(t)  +  Atf(t  +  — )  ( A2b) 

2 


A6 


This  scheme  is  accurate  to  second  order.  Applying  this  method 
(25): 


to  Eqrs 


A  7 


determine  { f) .  v  .  Also,  Eqn  (44a)  will  fix  {E}  .  Then,  Eqns  (44b) 
and  (43b)  calculate  (E)  ^  and  respectively.  The  whole  process 
can  now  be  repeated.  Note  that  since  E  and  f  are  determined  at  the  same 
time  intervals,  it  is  unimportant  which  order  they  are  computed. 

COMPUTATIONAL  IMPLEMENTATION.  The  one-dimensional  SGEMP  boundary 
layer  problem  has  now  been  reduced  to  a  set  of  algebraic  equations,  Eqns 
(43)  and  (44).  Now,  a  numerical  technique  must  be  picked  to  solve  these 
equations,  and  the  boundary  conditions  properly  applied. 

The  matrices  [A],  [A^],  and  tAvJj  are  NxN,  where  N  =  total  number 
of  nodal  points.  Since  there  are  N^  of  the  [A  »  the  problem  size  in 
terms  of  computer  memory  for  just  the  matrices  is  NxNx(N7z+2).  Thus,  if 
f>z=2Q  and  Nv=10,  the  rectangular  mesh  will  have  200  nodes,  with  the 
matrices  having  880,000  members!  Fortunately,  most  of  these  are  zero; 
that  is,  the  matrices  are  extremely  sparse.  In  fact,  it  can  be  shown 
that  for  the  rectangular  mesh, 


N , (A, A  )  =  (3N  -  2) ( 3N  -  2)  (45a) 

0  Z  Z  V 

N , (A  )  =  21N  N  -  14N  -  18N  +  12  (45b) 

where  N^  means  the  maximum  number  of  non-zero  members  for  the  specified 
matrices.  Also,  Eqn  (45b)  includes  all  of  the  NxN  matrices.  Thus, 
for  the  20  x  10  problem,  at  most  there  are  only  5376  non-zero  members. 

Because  of  this,  sparse  matrix  techniques  may  be  used  to  solve  the 
equations.  I  chose  a  fully  packed  scheme  which  stores  only  the  non-zero 
members  of  the  arrays.  Pointer  arrays  are  used  to  specify  which  row  and 
column  each  member  is  in.  The  algebraic  solution  is  done  using  a 
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Gauss-Seidel  iteration  technique  (Ref  35,  Sec.  20.8)  in  fully  packed 
form  (Ref  36).  Matrix  [A]  is  the  only  one  which  needs  to  be  inverted, 


and  is  diagonally  dominant  for  most  values  of  ^  .  Thus,  the 
Gauss-Seidel  method  is  guaranteed  to  converge. 

The  boundary  conditions  are  implemented  by  modifying  the  load 
vector,  {b}  ,  for  the  form  [A]  {x}  =  { b)  ,  as  well  as  the  appropriate 
row  and  column  of  [A]  (Ref  18,  Chapter  20).  This  reduces  the  number  of 
unknowns  in  {x}  by  the  number  of  boundary  values  needed  to  be  specified. 
Using  Eqn  (16)  and  (17),  the  following  boundary  condition  on  f  specifies 
the  source  of  electrons: 


f(0,v,t)  =  EJ.EX2)  T(t)  [sec/m4] 

1  2w,  J 


wi 


(46) 


v>0 


This  is  valid  for  an  exponential  energy  source. 

The  next  chapter  presents  the  results  obtained  from  the  solution  of 
the  equations  described  in  this  chapter  for  several  conditions. 


IV.  Results 


In  this  chapter,  the  numerical  results  obtained  by  the  application 
of  the  techniques  presented  in  Chapters  II  and  III  are  given.  These 
results  are  organized  in  order  of  increased  computational  difficulty. 
First,  sample  data  are  given  for  pure  advection,  or  wave  propagation. 
Then,  both  SGEMP  equations  are  solved  for  the  linear  case.  The  physical 
implications  are  analyzed  and  compared  with  the  computer  generated 
results.  Next,  the  nonlinear  set  of  equations  are  solved,  first  with  a 
linear  time  history,  and  then  with  an  exponential  time  dependence.  This 
is  compared  with  "particle  pushing"  techniques  and  quasi-static  theoret¬ 
ical  analysis.  Finally,  general  considerations  about  mesh  sensitivity, 
particle  conservation  and  convergence  criteria  shall  be  analyzed. 

PROPAGATION  STUDIES.  The  simplest  form  that  Eqn  (5)  can  have  and 
still  carry  some  physical  meaning  is  obtained  by  dropping  the  nonlinear 
term  to  give: 

|l(z,v,t)  +  v§I(z,v,t)  =  0  (47) 

3 1  3  z 

This  equation  describes  pure  advection,  or  propagation,  of  the  unknown 
function,  f(z,v,t),  at  velocity,  v.  It  is  still  a  first  order,  hyper¬ 
bolic,  partial  differential  equation,  but  is  now  linear.  Given  initial 
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conditions  that  prescribe  the  shape  of  the  distribution,  say  f  (z,v)  = 

o 

f(z,v,0),  then  the  exact  solution  to  Eqn  (47)  is  that  function  which 
translates  the  shape  of  f  at  the  velocity,  v,  along  the  z  axis. 
Therefore,  the  solution  to  Eqn  (47)  can  be  represented  by: 


f(z,v,t)  =  fQ(z  -  vt , v )  (48) 


This  problem,  and  numerical  methods  to  solve  it,  has  been  the 
subject  of  several  studies  (Ref  19,  Chapter  19  and  Ref  27).  In  particu¬ 
lar,  a  very  stringent  propagation  test  can  be  done  when  Eqn  (47)  is 
solved  using  discontinuous  initial  conditions  in  z.  A  square  wave  pulse 
meets  these  conditions. 

The  primary  purpose  in  doing  such  a  test  is  to  determine  how 
accurately  the  finite  element  technique  treats  a  propagating  signal. 
Since  oscillations  are  a  well-known  result  of  numerical  solutions  to 
wave  equations,  this  test  also  helps  to  study  the  effect  of  various 
upwinding  parameters  on  these  oscillations.  Several  computer  runs  were 
made  to  evaluate  different  .  The  test  was  done  in  a  manner  similar 
to  Book  and  Boris  (Ref  27).  The  pulse  existed  only  at  one  velocity, 
with  vftt/Az  *  0.2.  The  background  density  was  0.5  and  the  pulse  height 
was  2.0  (both  quantities  measured  in  arbitrary  units).  All  the  la^l, 
i=l,...,4  were  set  to  the  same  magnitude,  with  ct^  and  0*3  positive  for 
positive  v^  and  negative  for  negative  v^  .  Both  a2  and  04  were  always 
negative.  The  results  of  this  study  are  shown  in  Fig.  7,  after  100 
time  steps.  They  are  superimposed  upon  the  initial  square  wave.  Note 
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Fig.  7  Effect  of  Upwinding  on  Square  Wave  Propagation 


Che  characteristic  oscillations  in  front  of,  and  trailing,  the  pulse. 
It  is  clear  that  the  effect  of  increasing  the  magnitude  of  a  from  0  to  1 
is  to  reduce  the  oscillations,  while  introducing  dispersion  in  the 
pulse. 

There  is,  therefore,  a  trade-off  between  oscillatory  behavior  and 

dispersion.  A  value  of  about  0.5  appears  to  be  a  logical  choice  for  a 

compromise.  However,  for  practical  considerations,  another  value  makes 

more  sense  for  most  problems.  Consider  Eqn  (40a).  When  |cx  ^ I  =  2/3,  and 

the  conditions  on  the  sign  of  the  a  ^  are  met  as  described  above, 

a12=  a13=  a3 1  =  a32=  34 1  =  ai,2=0*  Because  six  out  of  16  members  of  this 

element  matrix  are  zero,  there  is  a  considerable  reduction  in  the  number 

of  non-zero  members  of  [A].  The  savings  is  slightly  more  than  a  factor 

of  two  for  a  regular,  rectangular  mesh.  Since  this  is  the  coefficient 

matrix  for  the  system  of  of  algebraic  equations  which  determines 

f(z,v,t),  execution  time  is  halved  by  using  | ct ^ |  =2/3,  i=l,...,4.  This 

is  due  to  the  decrease  in  the  number  of  arithmetic  operations  required 

by  the  Gauss-Seidel  iteration  scheme.  Only  the  non-zero  members  of  [A] 

are  used  for  the  tightly-packed  scheme  used  in  FEMNEP.  For  the 

square-wave  test,  there  are  847  non-zero  members  of  [A]  when  |a  J  =2/3. 

But,  when  |ot^|  *  2/3,  there  are  1810.  This  is  the  maximum  possible  for 

the  rectangular  mesh,  N^(A),  as  calculated  by  Eqn  (45a),  with  Nz  =61  and 

N  =4,  used  for  this  problem, 
v 

However,  an  interesting  feature  of  using  the  value  of  la^l  =2/3  is 
that  the  Gauss-Seidel  iteration  technique  does  not  converge  for  the 
square-wave  discontinuous  pulse.  For  all  other  values  of  a  tried,  the 
method  converged  for  this  problem,  including  a  =0.650  and  a  =0.670. 
Also,  the  iteration  scheme  worked  well  for  all  values  of  a^,  including 


=2/3,  for  smoothly  varying  input  pulses.  The  only  exception  to  this  ob¬ 
servation  occurs  when  |  a|  -*•  1 ,  at  which  point  the  Gauss-Seidel  technique 
no  longer  converges.  The  last  section  in  this  chapter  on  "Convergence 
Criteria”  gives  a  comparison  of  the  computer  execution  time  for  |a , | 
=2/3,  and  |a,|=0.60. 

ELECTRIC  FIELD  FROM  THE  LINEAR  VLASOV  EQUATION.  The  next  test  that 

can  be  made  on  the  finite  element  technique  is  to  add  the  electric  field 

3  f 

to  the  solution.  The  nonlinear  term,  E(z,t)- — (z,v,t)  is  still  ignored. 

a  v 

N, 

In  the  finite  element  approximation  this  implies  that  £*"  [A  ]  E  (t){f(t)} 

J-l  v  J  J 

is  zero  in  Eqn  (25).  Of  course,  the  physical  implication  of  making  the 

equations  linear  is  to  prevent  the  electrons  from  interacting  with  each 

other.  That  is,  an  electric  field  is  built  up  by  the  removal  of 

negative  charge  from  the  surface,  but  the  electrons,  themselves,  are  not 

affected  by  it.  The  linear  equations  allow  for  the  separation  of  charge 

with  no  boundary  layer  build  up.  In  effect,  a  parallel  plate  capacitor 

of  infinite  cross  sectional  area  is  the  physical  model  described  by  the 

linear  equations.  If  Q  represents  the  charge  per  unit  area  which  has 

been  emitted  from  the  surface,  then,  after  all  charge  has  proceeded  past 

distance  z  ,  the  electric  field  in  rationalized  MXS  units  is: 
o 

E(z)  =  -2-  =  ,  for  all  z<z  (49) 

G°  Eo 


and  N  =  number  of  electrons  per  square  meter. 

P 

The  purpose  of  running  this  test  is  to  check  the  effect  of  boundary 
conditions,  normalization  procedures  and  the  integrations  of  Eqn  (41). 
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The  test  was  performed  by  injecting  5.0  x  10^  electrons/cm^  into  a  4x20 
point  v-z  mesh.  The  electrons  were  given  only  one  velocity,  with  the 
condition  that  vAt/Az=0.2.  The  ct’s  were  all  set  with  the  same  magnitude 
for  each  element.  The  oq  and  013  were  positive,  while  the  0.2  and  at,  were 
negative.  Figures  8  and  9  show  the  results  of  this  test  for  two  cases, 

1 « ^ 1=2/3  and  a^=0.  Both  the  electric  field  vs  z,  and  the  distribution 
function  vs  z  are  shown  at  various  times.  The  incident  pulse  time 
history  was  of  the  form: 


T(t) 


a0 

a-B" 


(~e 


-ott 


+  e  ) 


(50) 


with  a  =  2/shake 
8  =  6/shake 

where,  1  shake  =  10  nanosec. 


Using  Eqn  (49),  the  theoretical  value  for  the  electric  field,  if 
6 

all  5  x  10  particles  have  passed  z„,  is  E  =  905.0  volts/m.  Table  I 

0  th 

shows  E  vs  z  produced  by  the  finite  element  equations  at  t=20.2 
shakes  (202  time  steps).  From  the  table  and  the  graphs,  it  is  easy  to 
see  the  effect  of  using  some  upwinding,  versus  not  using  any,  on  the 
electric  fir  .  The  effect  on  the  distribution  function  is  even  more 
dramatic.  From  Fig.  9,  not  only  are  trailing  oscillations  severe  with¬ 
out  upwinding,  but  reflections  off  the  end  of  the  mesh  at  z  =  95  cm  are 
are  very  apparent.  This  is  quite  visible  in  Fig.  9(f)  and  9(g),  where 


55 


TABLE  I 


Electric  Field  for  Two  Values  of  Alpha  in  a  Linear  Problem 


z  (  cm) 

E(z)  (volts/m) 
a  =  0 

E(z)  (volts/m) 
a  =  2/3 

0 

904.4 

904.4 

5 

913.5 

904.4 

10 

898.3 

904.4 

15 

889.9 

904.4 

20 

922.8 

904.4 

25 

906.5 

905.8 

30 

899.  7 

905.6 

35 

897.5 

906.4 

40 

899.2 

906.7 

45 

931.5 

906.6 

E  =  905.0  volts/m 
th 


the  reflected  distribution  function  is  growing  in  magnitude,  while  with 


upwinding,  Fig.  8(f)and  8(g),  reflections  do  not  exist.  This  wiil  be 


particularly  important  when  electrons  are  allowed  to  return  to  the  sur¬ 


face  with  negative  velocity.  Reflections  in  the  distribution  function 


from  returning  electrons  would  be  disastrous  from  three  standpoints: 


first,  positive  reflections  in  the  distribution  function  would  rep¬ 


resent  particles  external  to  the  object  which  are  not  really  there; 


second,  growing  reflections  represent  an  instability  in  the  calcu¬ 


lation;  and  third,  negative  reflections  giving  rise  to  a  negative 


distribution  function  have  no  physical  meaning  and  can  cause  instabili¬ 


ties  in  nonlinear  calculations. 


LINEARLY  RISING  PULSE.  The  next  type  of  solution  presented  is  for 


the  case  of  a  linearly  rising  pulse,  with  an  exponential  (black- 


body-like)  energy  spectrum.  This  is  a  truly  nonlinear  problem,  so  that 


the  full  set  of  equations  must  be  used.  As  a  check  on  the  results,  the 
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Fig.  8  (continued) 
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Fig.  9  (continued) 


finite  element  solution  is  compared  first  with  steady-state  theory  (Ref 
5).  This  is  permissible  because  whenever  the  photon  flux  rises  linearly 
with  time,  the  electron  behavior  can  become  quasi-static  under  certain 
conditions.  Then,  the  time  dependent  term  in  the  Vlasov  Equation  can  be 
ignored.  If  the  energy  spectrum  is  assumed  to  be  exponential,  the 
remaining  equations  can  be  reduced  to  quadratures  and  a  semi-analytic 
form  for  the  solution  produced. 

The  time  dependence  for  this  case  is  given  by  the  ramp  function: 


T(t) 


bt  0  <  t  <  t 

m 

0  <  t  <  t 


b  =  rise  rate  (sec-2) 


Unit  normalization  of  this  function  yields,  b=2/t2  .  In  order  to  run 

the  same  sample  problem  that  is  carried  out  in  Section  11  of  Ref  5,  the 

flux,  <t>Q ,  must  be  set  to  1/3  cal/cm  ,  with  tffl  =1  snake.  The  material  of 

the  object  is  assumed  to  be  Aluminum,  so  the  material  yield,  Y,  is  given 
12 

as  2.57x10  electrons/cal ,  again  taken  from  Ref  5.  Also,  the  exponenti¬ 
ation  energy  for  this  material  is,  wj  =4.77  keV. 

Using  these  parameters,  the  analysis  given  by  Carron  in  his  report 
will  be  valid  for  Eqn  (51)  for  0. 1  <  t  <  1  shake.  During  this  time, 
most  of  the  electrons  will  be  in  steady-state  conditions.  However,  all 
the  electrons  do  not  reach  steady-state  at  the  same  time,  since  they 
have  a  distribution  of  energies.  Thus,  caution  must  be  used  in  any 
detailed  comparison.  One-dimensional  SGEMP  can  be  solved  exactly  when 
all  the  electrons  are  In  steady-state,  but  the  application  to  any  type 
of  time  dependent  problem  must  be  considered  only  an  estimate.  Nonethe- 
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calcula- 


less,  it  is  instructive  to  compare  this  theory  with  a  dynamic 
tion  like  FEMNEP,  for  it  is  the  closest  one  can  come  to  an  analytic 
solution  for  a  nonlinear  problem. 

Figure  10  displays  an  example  of  what  a  typical  (z,v)  mesh  looks 

like  for  the  linearly  rising  pulse  problems.  It  is  now  necessary  to 

allow  particles  to  return  to  the  surface,  which  means  that  negative 

velocities  must  be  included  in  the  mesh.  The  sign  of  a  must  always  be 

i 

chosen  in  the  direction  of  motion  of  particles  through  the  mesh  in  both 
the  v  and  z  variables.  These  directions  are  shown  in  Fig.  10. 
Equation  (51)  and  Eqn  (46)  were  used  to  input  the  source  of  electrons. 

Several  different  mesh  sizes  were  tried,  and  compared  with 
steady-state  theory.  The  results  are  given  in  Fig.  11,  superimposed  on 
the  theoretical  predictions.  For  all  these  calculations , | a , | =  2/3.  In 
order  to  make  the  comparisons  easier,  the  meshes  were  set  up  using 
uniform  spacing  in  both  z  and  v  for  this  case  (except  near  zero  velocity 
where  20  cm/sh  =  1  IceV  is  the  smallest  velocity  allowed). 

Figure  11  clearly  shows  that  as  the  mesh  size  is  decreased,  conver¬ 
gence  is  obtained.  The  convergence  is  not  to  the  steady  state  values, 
however.  But,  this  is  to  be  expected.  Mote  that  close-in,  out  to  about 
three  Debye  lengths  (where  a  Debye  length  is  approximai  •  0.27  cm, 

according  to  Carron)  FEMNEP's  results  agree  to  withi*-  },  (x\  .all  that 

the  boundary  layer  is  defined  as  a  few  Debye  lengths  thick).  Figure  11 
also  shows  that  at  10  Debye  lengths,  steady-state  theory  and  the 
calculated  results  differ  by  more  than  a  factor  of  two.  Physically, 
this  is  expected.  The  higher  energy  electrons  which  penetrate  the 
boundary  layer  have  not  reached  steady-state  conditions.  But,  most  of 
the  lower  energy  electrons  which  are  trapped  in  the  layer  have  come  to 
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FEMNEP  VS  STEADY  STATE  THEORY 
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Fig.  11  Finite  Element  Comparison  with  Steady  State 
Theory  for  Linearly  Rising  Pulse 
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equil  ibrium. 

The  next  comparision  of  a  linearly  rising  pulse  is  made  with  a 
scaled  "particle-pushing"  program  (Ref  7),  called  SCAL1D.  Figures  12 
and  13  demonstrate  the  electric  field  comparison  obtained  as  a  function 
of  distance  and  time.  Both  figures  show  results  at  early  time,  that  is, 
before  1  nsec.  The  problem  being  solved  is  exactly  the  same  one  used  in 
the  last  paragraph.  The  mesh  scheme  is  similar  to  Fig.  10.  However, 
at  these  early  times,  the  electrons  are  not  in  equilibrium,  and  the  full 
time  dependent  solution  is  being  tested.  These  data  show  that  FEMNEP 
compares  very  favorably  with  SCAL1D.  For  all  times  and  distances  shown, 
the  agreement  is  no  worse  than  25%.  In  many  cases,  it  is  much  better 
than  that.  This  result  further  supports  the  conclusion  reached  in  the 
last  section  which  suggested  that  at  t  -  0.124  shake  many  high  energy 
electrons  have  not  reached  steady-state.  This  causes  the  agreement 
between  FEMNEP  and  steady-state  theory  to  deteriorate  beyond  several 
Debye  lengths. 

EXPONENTIAL  TIME  HISTORY.  One  final  comparison  has  been  made  using 
the  finite  element  computer  program  in  a  calculation  with  an  exponential 
time  history,  Eqn  (18).  This  time,  the  results  are  compared  to  another 
"particle-pushing"  program  called  MAD1  (Ref  8).  The  major  difference 
between  this  code  and  SCAL1D  is  that  the  MAD1  code  inputs  the  energy 
spectrum  of  the  electrons  from  a  one-dimensional  emission  code,  rather 
than  treating  it  as  an  analytic  form.  Of  course,  none  of  the  unknown 
variables  are  scaled  in  '1AD1  as  they  are  in  SCAL1D. 

The  problem  parameters  used  for  this  comparison  are  as  follows:  a  5 
keV  blackbody  photon  spectrum,  normally  incident  on  Aluminum  with  a  time 
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FEMNEP  VS  SCAL1D 


Fig.  12  Finite  Element  Comparison  with  Scaled  Results, 
Electric  Fields  Vs  Distance 
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.  13  Finite  Element  Comparison  with  Scaled  Results, 
electric  Field  Vs  Time 


history  defined  by: 


cj  =  2.0/shake 
C£  =  2.5/shake 
c3  =  1.569/shake 
c4  =  4.0  shakes 

and,  the  total  X-ray  fluence,  <Ji  =  0.1  cal/cm2. 

For  the  finite  element  calculation,  the  same  source  parameters  were 

12 

used  as  in  the  previous  case;  that  is,  Y  =  2. 57  x  10  elect rons/cal  and 
wi  =  4.77  keV,  assuming  that  the  blackbody  radiation  is  represented  by 
an  exponential  energy  distribution  of  electrons. 

Using  the  identical  upwinding  parameters  as  before,  la. I  =2/3,  for 
all  i,  the  finite  element  calculation  was  performed  in  a  nonuniform 
grid.  This  was  necessary  so  that  distances  out  to  five  meters  could  be 
used,  while  still  resolving  the  Debye  length  close  to  the  surface. 
Table  II  summarizes  all  the  input  parameters  for  both  the  MAD1  and 
FEMNEP  runs.  It  is  interesting  to  note  that  the  ‘1AD1  calculation  was 
done  with  a  uniform  z-spacing  of  1.5  cm.  The  FEM  required  smaller 

spacing  in  close  (see  next  section  on  mesh  sensitivity  for  discussion  of 
z-spacing  requirements  for  FEMNEP). 

Figures  14  and  15  show  the  comparison  of  the  electric  field  versus 
time  on  the  surface,  and  electric  field  versus  distance  at  t=4.0  shakes. 
The  discrepancy  at  very  early  time  in  Fig.  14  is  due  to  the  fact  that 
MAD  1  must  start  its  calculation  at  t=-2  shakes  to  avoid  injecting  too 
many  particles  in  the  beginning.  The  agreement  over  most  of  the  z  and  v 
range  is  excellent. 

MESH  SENSITIVITY  STUDIES.  The  purposes  of  this  section  are 
twofold:  (1)  to  describe  the  important  restrictions  and  requirements  on 
mesh  spacing  and  time  step  sizes,  and  (2)  to  demonstrate  the  convergence 
of  the  FE!t  as  mesh  size  decreases.  Examples  are  given  which  show  that 
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TABLE  II 


Input  Parameters  for  MADi  VS  FEMNEP  Comparisons 


MAD1 


FEMNEP 


CLASSIFICATION 

NUMERICAL  METHOD 

SOURCE  ENERGY 
SPECTRUM 

SOURCE  ANGULAR 
DISTRIBUTION 

NO.  OF  EMISSION 
ENERGY  BINS 

NO.  OF  EMISSION 
ANGULAR  BINS 

EON  USED  FOR 
E-FIELD  CALCULATION 

E-SPACING 

AZ 

TOTAL  Z  CELLS 
TIME  ALGORITM 

At 

NO.  OF  TIME  STEPS 
TO  6.0  SHAKES 


Particle  Pusher 

Finite  Differences 

From  QUICKIE2  Code 
(Ref  31) 

Assumed  cos  0 

20 

9 

Gauss'  Law 
Uniform 
i .  5  cm 

400 

2nd  Order 
centered 

finite  difference 
0.01  shake 

800 

(calculation  starts 
at  t  =  -2.0  shakes) 


Vlasov  Eqn  Solver 
Finite  Elements 
Assumed  Exponential 

Assumed  cos  9 

5 

Analytic 

Ampere's  Law 

Non-uniform 

0. 18  cm  minimum 
30  cm  maximum 

44 

2nd  Order 
Lax-Wendrof f 

0.001  shake 

6000 
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Fig.  14  Finite  Element  Comparison  with  MAD1  Program, 
Electric  Field  Vs  Time 
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Fig.  15  Finite  Element  Comparison  with  MAD1  Program, 
Electric  Field  Vs  Distance 


certain  recommendations  must  be  followed  to  preserve  the  integrity  of 
the  finite  element  solutions.  And,  results  will  be  presented  indicating 
the  finite  element  solutions  converge,  for  decreasing  mesh  sizes,  to  the 
same  data  that  a  particle  code  produces. 

The  first  area  of  concern  are  the  Courant-Friedrichs-Lewy 
conditions  of  Eqns  (10)  and  (11).  They  impose  restrictions  on  the 
relative  sizes  of  the  Az  and  At  spacing.  Figures  16  and  17  show  what 
occurs  to  the  electric  field  and  distribution  function  when  At  is 
increased  from  0.0001  shake  to  0.002  shake  for  a  fixed  Av  and  Az  mesh  of 
30  cm/sh  and  0.1  cm,  respectively.  This  is  a  variation  of  0.09  to  1.8 
for  the  Courant  factor,  vAt/Az  when  v=90  cm/sh  (the  largest  velocity 
in  the  mesh  which  particles  were  allowed  to  travel).  The  results  show 
vividly  that  the  Courant  condition  of  Eqn  (14)  must  be  adhered  to 
strictly  to  preserve  stability.  The  cross-over  from  stability  to  insta¬ 
bility  occurs  at  At=0. 001111  shake,  or  vAt/Az=l.  The  results  also  show 
that  once  vAt/Az  becomes  less  than  one,  stability  is  quickly  achieved, 
and  no  matter  how  small  At  is  taken  to  be,  there  is  no  noticeable  change 
in  the  electric  field  profiles.  There  is,  however,  an  improvement  in 
particle  conservation,  which  will  be  taken  up  in  the  next  section.  An 
indication  that  there  is  a  net  effect  by  making  At  small  enough  can  be 
seen  in  Fig.  17.  Oscillations  of  the  distribution  function  are  seen  at 
the  highest  velocity  profile  (v«90  cm/sh)  for  large  At. 

A  comment  should  be  made  at  this  point  about  the  other  Courant 
condition,  Eqn  (11).  This  condition  can  be  re-written  in  terms  of  the 
electric  field,  with  the  help  of  Eqn  (9): 

AtSshl"^'  >  1.76  10_  ?E[  volts/m]  (52) 


c n  m 
co 

^  > 
T? 

3  U 
4-»  O 
CO 


i 


77 


d(V/m)  Vs  Distance(cm) 


(e) 


16  (continued) 
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Fig.  17  (continued) 
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The  above  equation  will  be  satisfied  for  all  practical  situations.  The 
Av  spacing  will  not  be  much  smaller  than  10  cm/sh=0.29  keV,  and  at  is 
typically  about  0.01  shake  or  less.  Therefore,  the  smalxest  that  the 
left-hand-side  of  Eqn  (52)  would  ever  get  is  about  1000  cm/sh  .  In 
order  for  Eqn  (52)  not  to  hold  true^  E  must  be  larger  than  1.0  x  101() 

volts/m - very  large  indeed.  Also,  very  conservative  figures  have  been 

used  to  determine  Av/At.  For  example,  the  calculations  presented  in 
this  dissertation  used  a  Av/At  closer  to  20,000  cm/sh  .which  pushes  the 
critical  value  of  E  even  higher.  Consequently,  this  Courant  condition 
is  easily  satisfied  for  all  practical  one  dimensional  SGEMP  problems. 

The  importance  of  choosing  a  sufficiently  small  Az  spacing  was 
discussed  in  Chapter  II.  The  physical  nature  of  the  Debye  length 
requires  Az  to  be  able  to  resolve  it.  One  way  of  experimentally 
measuring  the  Debye  length  in  a  calculation  is  to  take  it  as  the 
distance  the  electric  field  falls  by  a  factor  of  2/e.  This  comes  from 
the  screening  factor  of  exp(-z/AQ)  in  the  static  equations  for  the 
electric  potential  (Ref  28).  From  Fig.  II,  E  falls  from  14.5  x  I05 
volts/m  to  10.5  x  io5  volts/m  (factor  of  2/e)  in  about  0.15  cm.  If  the 
z-spacing  exceeds  0.2  cm,  the  results  should  be  suspect.  The  same 
problem  as  shown  in  Fig.  11  was  run  with  a  z-spacing  of  0.4  cm.  These 
data  are  presented  in  Fig.  18,  along  with  the  steady-state  predictions, 
and  the  Az=0.20  run  for  comparison.  Not  only  is  the  curve  for  Az=0.40 
cm  less  accurate,  as  might  be  expected,  but  the  shape  of  the  curve  shows 
oscillations.  In  fact,  at  distances  beyond  3  cm  ,  E  becomes  negative 
and  oscillating. 

The  restrictions  on  the  Av  and  Az  mesh  sizes  mentioned  above  are 
the  minimum  requirements  in  relation  to  At  step  size.  However,  there 
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Fig.  18  Finite  Element  Calculation  Showing  the 
Effect  of  Az  Spacing  Large 
Compared  to  Debye  Length 
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TABLE  III 


Electric  Field  at  Various  Distances 
Under  Mesh  Refinement  at  t  =  0.08  shake 


RUN# 

Az 

(cm) 

Av 

(cm/sh) 

E(0) 

(MV/m) 

E(. 257cm) 

(MV/m) 

E( . 857cm) 

(MV/m) 

1 

.05* 

30 

.9670 

.6794 

.3438 

2 

.05* 

20 

.9304 

.5694 

.2516 

3 

.05* 

15 

.8998 

.5105 

.2136 

4 

.05* 

10 

.8573 

.4605 

.1932 

5 

.025 

10 

.8526 

.4600 

- 

6 

.05 

5 

.8093 

.4355 

- 

7 

.05 

1* 

.7950 

.4477 

- 

SCAL1D 

— 

— 

.8455 

.4585 

.1876 

(MV/m  =  Megavolt/meter) 

*  non-uniform  mesh,  number  represents  the 
smallest  mesh  interval 


are  other  considerations  equally  as  important  with  regard  to  the  conver¬ 
gence  and  accuracy  of  the  method.  In  order  to  demonstrate  that  the  FEM 
does  produce  results  which  converge  to  other  known  data,  a  series  of 
runs  were  made  in  which  the  mesh  spacing  was  systematically  reduced. 
The  data  were  then  compared  with  the  most  accurate  and  detailed  particle 
code  results  available  for  one-dimensional  SGEMP,  the  SCAL1D  program. 
Table  III  summarizes  the  results  of  this  study. 
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The  data  in  Table  III  are  analyzed  the  following  way:  as  the  v 
mesh  size  decreases,  the  electric  fields  approach  converging  values.  At 
the  last  refinement,  the  difference  in  the  calculated  fields  does  not 
exceed  6%,  and  over  much  of  its  range,  it  is  better  than  4%.  Also,  for 
the  two  values  for  which  the  FEtlNEP  data  can  be  compared  with  SCAL1D 
results  directly,  without  interpolation,  they  lie  within  6.5%.  The 
crossover  of  the  bottom  two  curves  is  due  to  the  non-uniform  spacing  in 
the  z-direction  for  this  run.  A  graphical  representation  of  these 
results  appear  in  Fig.  19.  When  this  data  is  coupled  with  the  z-mesh 
refinement  studies  of  Fig.  11,  convergence  of  this  finite  element 
approach  to  the  one-dimensional  SGEHP  boundary  layer  is  shown.  Also, 
the  comparison  of  FEMNEP  data  with  the  particle  code  SCAL1D  demonstrates 
the  validity  of  the  technique. 

As  discussed  above,  the  Courant  condition  for  the  v-mesh  does  not 
apppear  to  play  a  significant  role  in  the  choosing  of  a  Av  since  it  is 
normally  satisfied  for  practical  SGEHP  problems.  However,  Fig.  20 
shows  that  the  v-step  can  not  be  picked  arbitrarily  small,  without 
regard  to  the  z-step.  The  build-up  of  the  knee  in  the  curve  for  Az  = 
.05  cm  is  not  normal.  It  is  a  result  of  choosing  too  small  a  Av  for 
this  Az.  Note  that  the  knee  is  not  present  for  the  Az  =  .05  cm,  Av  =  20 
cm/sh  curve.  (This  curve  diverges  from  the  other  two  curves  beyond  z  = 
.15  cm  because  a  different  boundary  condition  was  used  at  z  =  .5  cm.) 
The  problem  is  corrected  by  making  Az  smaller.  This  coupling  effect 
between  Az  and  Av  appears  to  get  more  severe  as  the  number  of  time 
iterations  increases,  and  impacts  the  flexibility  of  grid  make-up  for 
problems  with  long  run  times. 

PARTICLE  CONSERVATION.  The  computer  program,  FEMNEP,  keeps  track 
of  the  number  of  electrons  emitted  at  each  time  step,  and  the  number 
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•  =  FEMNEP  DATA 

<S>  =  SCAL1D  DATA 


0  0.1  0.2  0.3  0.4  0.5 

z  (cm) 


Fig.  19  Electric  Field  Vs  Distance  at  t  =  0.08  sh, 
Showing  Convergence  to  SCAL1D  Results 


z  (cm) 


Fig.  20  Electric  Field  Vs  Distance  at  t  =  0.16  sh, 
Showing  Effect  of  Small  A v  Spacing 
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leaving  (or  entering)  the  mesh  from  the  free  boundaries  (see  Fig.  1). 
The  total  number  of  electrons  in  the  mesh  can  also  be  calculated.  From 
this  information,  particle  conservation  is  checked.  The  parameters 
needed  are: 


MP  =  total  particles  in  mesh  per  cm2 

MR  =  particles  returned  to  surface  per  cm2 

NL  =  particles  leaving/entering  mesh  at  z  per  cm2 

max 

MB  =  particles  leaving/entering  mesh  at  -v  per  cm2 

max 

PE  =  particles  emitted  into  mesh  per  cm2  . 


Then,  PB  =  particle  balance  =  PE+NR-NL+NB.  For  perfect  conservation,  PB 
would  equal  MP  at  all  times  in  the  calculation. 

These  parameters  are  determined  from  the  following  equations: 


MP  = 


z  +v 
max  max 

/  Je(z-v- 


t)  dvdz  (53) 


0  -v 


max 
t  +v 


NR 


a  m  max 
I  /vf (0 , v, t ' 


)  dvdt'  (54) 


NL 


//■ 

0  0 


max 

+v 

max 

vf(z  ,v,t')  dvdt'  (55) 
max 
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NB 


t  z 


-?/  /e<2' 

0  0 

t  V 

.  .  max 

PE  =  /  /vf(0,v,t' 


t,)f(z‘~Vmax‘t,)  dzdt' 


)  dvdt '  (57) 


(56) 


The  first  equation  is  the  definition  of  the  distribution  function,  and 
the  other  equations  are  expressions  for  the  flow  rate  of  particles 
across  a  boundary,  integrated  over  time. 

Equation  (53)  can  be  expressed  as  a  simple  sum  over  the  elements  in 
the  rectangular  mesh: 


e  4 

NP  =J\(e)b(eQ[Ve)(t:)  1  (58) 

e=»l  i=l 

where  f  (t),  1*1,... ,4,  are  the  four  nodal  values  of  the  distribution 
i 

function  for  element  (e)  at  time,  t.  The  other  equations  are  in  the 
same  form  as  the  integral  in  Eqn  (6).  Therefore,  Eqns  (54)  thru  (57) 
are  calculated  using  the  procedure  described  by  Eqn  (41).  For  example, 
Eqn  (57)  becomes, 


PE  -  N  (At  )fS  \v  -v  )[2v  f  (0 , v  ,t  )  +  v  f (0,v  ,t  ) 

1 Z-J  k+1  k  k  k  i  k  k+1  t 


i-1  k-k 


+  v  f ( 0 , v  ,t  )  +  2v  f ( 0 , v  ,t  )]  }  (59) 

k+1  k  i  k+1  k+1  i 
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where, 


N  =  number  of  time  steps  out  to  t 
t 

kQ  =  velocity  subscript  assigned  to  v  =  0 

k  =  velocity  subscript  assigned  to  v  =  v 
m  max 

The  other  equations  are  computed  in  a  like  manner. 

Figure  21  shows  these  parameters  for  a  typical  linearly  rising 
pulse.  These  data  correspond  to  the  same  output  shown  in  Fig.  13.  In 
the  graph,  PE  =  emitted  curve,  NR  =  returned  curve,  NP  =  total  curve, 

and  NB  +  NL  =  lost  curve.  This  means  that  the  curve  labeled  "lost" 

could  actually  represent  a  gain.  However,  its  real  meaning  is  the  total 
number  of  electrons  which  either  leave  the  mesh  or  enter  it  from  a  free 
boundary.  This  figure  is  an  excellent  presentation  of  the  overall 

behavior  of  the  electrons  as  a  function  of  time.  It  shows  that  at  early 

time,  all  the  particles  emitted  stay  in  the  mesh.  Eventually,  the 
electric  field  returns  so  many  electrons  to  the  surface  that  the  total 
number  of  them  levels  off.  It  also  shows  the  number  of  particles  "lost" 
to  the  calculation  is  insignificant  with  respect  to  all  other  particle 
paramters. 

Figure  22  shows  the  same  information  for  the  exponentially  rising 
pulse  corresponding  to  the  electric  fields  of  Figs.  14  and  15.  In  this 
case,  the  curves  clearly  show  that  equilibrium  between  the  emitted  and 
returned  electrons  only  occurs  near  the  peak  of  the  pulse  at  t*4.0 
shakes.  Before  this  time,  and  shortly  after  5  shakes,  the  total 

particles  in  the  mesh  is  a  strong  function  of  time.  Also,  it  is  clear 

that  lie  number  of  electrons  lost  to  the  mesh  plays  no  role  whatsoever. 

The  next  two  figures  present  examples  of  the  efficiency  of  particle 
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ELECTRONS/CM 


TIME  (SHAKES) 


Fig.  22  Particle  Conservation  Parameters  for 
Exponentially  Rising  Pulse 
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conservation  for  several  finite  element  calculations.  Particle  conser¬ 
vation,  PC,  is  defined  as: 


PC 


100  | NP-PB | 
NP 


(60) 


The  smaller  this  number  is,  the  better  the  conservation  of  electrons  is 
preserved.  Figure  23  shows  PC  as  a  function  of  time  for  FEMNEP  in  the 
MAD1  comparisons.  For  this  particular  calculation,  electrons  are 
conserved  to  an  accuracy  of  better  than  90%.  It  is  interesting  to  note 
the  tracking  of  this  parameter  with  the  source  time  history. 

The  final  graph  in  this  chapter,  Fig.  24,  displays  the  effect  of 
decreasing  the  time  step  on  particle  conservation.  The  problem  used  in 
this  example  is  the  linear  time  history,  and  results  of  Figs.  16  and 
17.  Over  most  of  the  time  range  shown,  the  electron  conservation  is 
better  by  a  factor  of  2  to  3  for  the  smaller  t.  The  anomolous  dips  in 
both  curves  are  a  curious  aspect  of  this  procedure  which  also  shows  up 
in  the  exponential  case  of  Fig.  24. 

CONVERGENCE  CRITERIA.  The  Gauss-Seidel  algorithm  is  used  as  the 
algebraic  equation  solver.  An  absolute  convergence  test  is  applied  to 
determine  whether  another  iteration  is  required.  That  is,  iterations 
will  continue  until, 


(v+1)  (v) 

xi  '  xi 


<  e 


(61) 
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TIME  (SHAKES) 


Fig.  23  Particle  Conservation  for  the  flADl  Comparison 


where , 


the  ith  unknown  after  the  vth  iteration 


(v) 


£  =  absolute  convergence  criteria 


The  initial  guess  for  x^^is  taken  to  be  the  last  calculated  values  for 
at  the  previous  time  step,  (zero  for  the  first  guess).  In  all  of  the 
results  presented  thus  far,  £  =  0.001,  This  is  a  very  stringent  test 
because  the  unknown  function,  f(z,v,t),  is  of  the  order  of  1.0  x  105  for 
these  problems.  Table  IV  shows  the  effect  of  relaxing  the  convergence 
test  on  the  same  problem  solution  for  three  different  values  of  £  . 

Since  program  execution  time  listed  in  the  table  includes  entire 
program  execution,  the  percent  reduction  in  running  time  will  increase 
slightly  for  longer  running  jobs  (see  below).  From  the  data  in  Table 
IV,  it  is  obvious  that  a  much  less  severe  convergence  test  than  0.001  is 
perfectly  acceptable.  It  is  also  clear  that  this  factor  has  an 
important  impact  on  computer  costs.  Since  no  special  attempts  were  made 
to  optimize  execution  time  for  the  finite  element  method,  comparisons 
with  other  techniques  to  determine  the  "fastest"  method  are  not 
possible. 

Note  that  when  |aj  ^  2/3,  the  execution  time  increases  by  55%  com¬ 
pared  to  the  case  when  |a^|  =  2/3,  for  the  same  problem.  At  the  same 
time,  the  number  of  iterations  remain  about  the  same. 

The  large  e  was  used  on  several  different  cases  in  order  to  deter¬ 
mine  the  reduction  of  computational  time  which  can  be  realized  without 
affecting  the  electric  field  results  significantly  (2nd  or  3rd  signifi¬ 
cant  digit).  Table  V  shows  these  comparisons  for  different  size  jobs. 
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TABLE  IV 


Convergence  Parameters  for  a  Linearly-Rising 
Pulse  Problem  at  t  =  0.08  shake 


e 

EXECUTION 

TIME(SEC) 

f (z  ,v 
( sh?cm) 

)  E(0) 
MV/M 

TOTAL  ELECTRONS 
IN  MESH  (cm-2) 

ITERA¬ 

TIONS 

ft  N0N- 
ZEROS> [A] 

'«i 

0.001 

137.3 

2496 

.9635 

5.32643  x 1C9 

30 

728 

2/3 

100 

113.3 

2496 

.9635 

5.32643  x 109 

22 

728 

2/3 

100000 

91.1 

2465 

.9650 

5.32570  x io9 

13 

728 

2/3 

100000 

141.4 

2903 

.9646 

5.  33002  x  io9 

11 

1540 

.60 

zo 

=7.1  cm,  vQ  =  90  cm/sh 

TABLE  V 

FEMNEP  Execution  Time  and  Memory  Requirements 
For  Several  Different  Cases  (a=2/3) 


PROBLEM 

ft  OF 

EXECUTION 

%  REDUCTION 

MEMORY 

At  ft  OF 

TYPE 

£ 

NODES 

TIME(SEC) 

IN  EXECUTION 
TIME 

FOR  MATRIX 
STORAGE 

TIME 

(sh)  STEPS 

Steady-  .001 
State 

Comparison 

Fig.  11, 

128 

8.376 

8656 

.001 

125 

Av=30.  100000 

128 

6. 154 

27 

8656 

.001 

125 

Steady-  .001 
State 

492 

111. 183 

— 

35,054 

.0005 

250 

Comparsion 

Fig.  11, 

A v=20  100000 

492 

76.424 

31 

35,054 

.0005 

250 

SCAL1D  .001 

192 

137.313 

_ _ 

13,184 

.0002 

1204 

Comparison 

100000  192 


91.108 


34 


13,184  .0002  1204 

13,184  .0002  1204 


Memory  requirements  listed  in  the  table  are  for  [A],  [  Az  ]  ,  and  the 
NTz  [Av]  matrices,  and  their  row  and  column  pointer  arrays.  It  does  not 
include  other  miscellaneous  arrays  or  the  remaining  program  storage.  If 
lal#2/3  for  these  runs,  the  storage  for  matrix  [A]  would  increase  by  a 
factor  of  two,  and  execution  time  would  suffer,  as  shown  in  Table  IV. 
When  I ct  | =2/3  is  used,  there  is  a  slight  increase  in  the  percent  savings 
in  execution  time,  for  £  =  100000,  as  the  number  of  time  increments  get 
larger. 

The  FEMNEP  comparison  with  the  flADl  code,  Table  II,  took  2607 
seconds  (43.5  minutes)  to  execute  with  £  =0.001.  However,  since  6000 
time  steps  were  done  in  this  run,  a  factor  of  34%  savings  can  conserva- 
tivly  be  estimated  for  £  =  100000.  This  would  reduce  running  time  to 
about  720  seconds  (29  minutes).  Memory  requirement  for  the  matrices  in 
this  run  was  37,655  decimal  words,  for  I® I  =2/3,  and  528  nodes. 
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V.  Summary  and  Conclusions 


The  overall  goal  of  this  study  was  to  investigate  the  use  of  the 
finite  element  method  for  the  solution  to  the  one -dimensional  SGEMP 
boundary  layer  problem.  This  feature  of  SGEMP  plays  a  significant  role 
in  the  final  amount  of  current  which  flows  on  the  surface  of  the  space 
object.  Because  of  this  important  effect,  the  boundary  layer  has  been 
treated  using  many  different  techniques.  Numerical  approaches  are  the 
most  prevalent  methods  currently  being  used. 

SUMMARY.  The  one-dimensional  boundary  layer  was  described  by  the 
Boltzmann  transport  equation  for  electrons  using  an  analytic  source 
function.  This  equation  is  a  time-dependent,  nonlinear,  first-order, 
partial  integro-dif ferential  equation  in  three  independent  variables:  a 
space  variable,  z;  a  velocity  variable,  v;  and  time,  t.  Since  this  is 
an  initial  value  problem,  a  time  marching  scheme  can  be  used  effective¬ 
ly.  The  temporal  behavior  of  the  electrons  is  treated  independent  of 
the  sptce  and  velocity  variables. 

The  finite  element  equations  were  developed  for  the  (z,v)  variables 
using  a  regular  rectangular  mesh  and  linear  approximations.  The  Method 
of  Weighted  Residuals  was  used  to  derive  the  integral  relationship.  The 
choice  of  the  weight  functions  was  dictated  by  the  dominance  of  the 
advective  terms  in  the  transport  equation,  and  its  nonlinear  nature. 
Weight  functions  which  were  identical  to  the  approximating  trial 
functions  did  not  work.  Therefore,  weights  were  picked  which  depend  on 
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This  "up- 


the  electrons'  direction  of  travel  through  the  (z,v)  grid, 
winding"  technique  creates  a  new  parameter  which  determines  the  amount 
of  bias  to  be  applied  to  each  element.  For  these  equations,  a  minimum 
•number  of  non-zeroes  occur  in  the  coefficient  matrix  for  the  resulting 
set  of  algebraic  equations  if  a  specified,  constant  value  of  2/3  is  used 
for  the  magnitude  of  the  upwinding  parameter.  This  value  also  provides 
a  reasonable  compromise  between  accuracy  of  pulse  propagation  and 
dispersion  of  the  pulse. 

With  the  upwinding  scheme,  special  care  was  required  for  the  time 
marching  algorithm.  Second  order,  centered,  finite  difference 

techniques  were  not  successful.  However,  a  two-step  Lax-Wendrof f  method 
was  found  compatible  with  the  fin  ;e  element  portion  of  the  solution. 

Boundary  conditions  on  the  transport  equation  in  phase  space  were 
applied  in  a  manner  similar  to  neutral  particle  transport  techniques  in 
real  space.  With  the  upwinding  set  up  in  the  direction  of  electron 
motion  through  the  grid,  the  distribution  of  electrons  does  not  get  nu¬ 
merically  reflected  off  any  of  the  sides  of  the  mesh. 

The  algebraic  equations  were  solved  twice  every  time  step  using  an 
iterative  technique.  A  tightly  packed  scheme  was  chosen  which  only 
stored  non-zero  members  of  the  arrays.  This  method  converged  for  every 
choice  of  upwinding  parameter,  cx,  used,  with  the  exception  of  I ct  1  =  1. 

These  methods  and  procedures  were  used  to  generate  a  new  set  of 
results  for  the  electric  fields  developed  near  a  plane  Aluminum  surface 
exposed  to  X-rays.  These  finite  element  solutions  were  then  compared  to 
analytic  estimates  and  particle  simulation  data.  The  latter  comparisons 
involved  two  different  sets  of  data:  scaled  curves  based  upon  a  linearly 
rising  source,  and  a  calculation  for  an  exponentially  rising  and  falling 
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source.  Although  completely  different  approaches  to  the  same  problem 
were  used  for  both  cases,  very  good  quantitative  comparisons  were 
obtained. 

Finally,  many  computations  were  made  to  determine  the  behavior  and 
stablilitv  of  the  finite  element  results  as  a  function  of  the  increments 
used  for  all  three  independent  variables,  z,  v,  and  t.  Both  physical 
and  numerical  limitations  were  investigated. 

CONCLUSIONS.  The  finite  element  equations  have  been  developed  for 
application  to  the  SGEMP  boundary  layer  problem.  SGEMP  is  a  member  of  a 
class  of  charged  particle  transport  phenomena  which  are  described  by 
time-dependent,  nonlinear,  partial  integro-dif f erential  equations.  The 
results  of  this  study  indicate  that: 

(1)  Linear  trial  functions  can  be  used  successfully  with  a  weighted 
residuals  approach  if  the  weights  are  dependent  on  the  direction  of 
particle  travel  through  the  mesh. 

(2)  When  linear  weight  and  trial  functions  are  used,  a  value  of  2/3 
for  the  upwinding  parameter  reduces  the  storage  requirements  for  the  co¬ 
efficient  matrix,  [A],  when  iterative  solutions  are  sought.  It  also 
decreases  the  computational  time  needed  for  iterative  methods. 

(3)  The  FEM  can  be  implemented  for  the  transport  equation  and  the 
coupled  equation  for  the  electric  field  in  the  same  calculational  mesh 
when  linear  trial  and  weight  functions  are  chosen. 

(4)  Free  surface  boundary  conditions  can  be  use-4  in  the  computa¬ 
tional  mesh,  even  at  the  interface  where  electrons  are  returning  to  the 
surface  of  the  object. 

(5)  Courant  stability  conditions  are  sufficient  to  ensure  that  the 
solutions  to  the  equations  do  not  grow  exponentially.  This  is  true  for 
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both  the  z  and  v  variables.  However,  v  cannot  be  made  indiscriminately 
small  without  reducing  z  as  well. 

(6)  The  z  increment  must  be  chosen  small  enough  to  resolve  the 
plasma  Debye  length  with  several  mesh  steps. 

The  comparisons  of  the  finite  element  calculations  with  other 
methods  show  that  this  technique  can  be  used  to  analyze  the  difficult 
nonlinear  SGEMP  boundary  layer  with  success.  They  also  provide  new 

results  which  support  particle  simulation  methods  and  theory. 
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APPENDIX  A 


Derivation  of  the  One-Dimensional  Equations 

In  a  fully  three  dimensional  problem,  the  equations  governing  the 
notion  of  the  electron  plasma  are  given  by  Eqns  (1),  (2),  (3)  and  (4). 
When  a  reduction  to  one  spatial  dimension  and  one  velocity  direction  is 
desired,  the  assumption  is  made  that  variations  of  the  dependent 
variables  are  allowed  only  in  these  directions.  That  is, 


f(x,v,t)  = 

f(z,v  ,t) 
z 

>  -►  + 

a(x, v, t )  = 

a(z,v  ,t) 
z 

E(x,t)  = 

E(z,t) 

B(x,t) 

B(z,t) 

Figure  25  depicts  the  geometry  for  the  one-dimensional  problem.  Under 
these  conditions,  Eqn  (I)  immediately  becomes: 


3f 


3f 


3  f  , 


3T(z*Vc)  +  vz37(z’VC)  +  az37(z*vz’t)  =  S(z’V0 


(A-I) 


Thus,  only  a  (z,v  ,t)  plays 
z  z 


role,  and  it  is  now  reduced  to: 


a  (z,v  ,  t)  -  -  — E  (z,t)  (A-2) 

z  z  m  z 
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since  a  cannot  be  a  function  of  v  or  v  •  Also,  the  two  Maxwell's 
z  x  y 

Equations,  Eons  (3c)  and  (3d)  become: 


£?<z.t)  -f*(z,t)  (A-3a) 

p*(z,t)  =  -  ||y(z,t)  (A-3b) 


iiz 

at 


(z,t) 


=  o 


(A-3c) 


And, 


-  - 


3  E 

uoJx(z’C)  + 


(A-4a) 


H^(z,t)  =  u0Jy(z,t)  +  Uoeo|fy(z,t) 


(A-4b) 


J  z(  z  ,  t ) 


3E 
Eo3 1 


:(z,t) 


(A-4c) 


All  but  one  of  these  equations  can  be  eliminated  from  further  con¬ 
sideration.  For  example,  the  existence  of  E  (z,t)  will  create  an  accel¬ 
eration  in  the  x-directlon  from  Eqn  (2),  unless. 


Ex(z,t)  =  vzBy(z,t)  ”  vyBz(z,t)  (A-5) 
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But,  if  this  were  the  case,  Eqn  (A-5)  can  be  twice  differentiated  with 
respect  to  time,  to  give: 


(A-6) 


where  Eqns  (A-3b)  and  (A-3c)  have  been  used.  However,  when  Eqn  (A-4a) 
is  differentiated  with  respect  to  t,  and  Eqn  (A-3b)  is  substituted  in, 
the  result  is. 


x(z,t) 


1  32E 


x(z,t) 


=  0 


(A-7) 


since  Uoeo  =  1/c  • 

Comparing  Eqns  (A-6)  and  (A-7),  the  only  way  both  can  hold  true  is 
if  v^=  c.  Since  this  is  not  possible,  the  only  conclusion  that  can  be 
reached  is  that  Ex(z,t)=0.  In  a  similar  manner,  it  is  an  easy  task  to 
show  chat  E  (z,t)=0. 

Once  and  Ey  have  been  identified  as  zero,  Eqns  (A-3)  require 
that  3(z,t)  =  0,  because  B(z,0)  =  0.  The  only  remaining  Maxwell's 
Equation  is  Eqn  (A-4c),  Ampere's  Law.  This  is  one  of  the  two  fundamen¬ 
tal  equations  used  in  this  dissertation,  the  other  being  Vlasov's 
Equation,  Eqn  (A-l)  with  the  source  specified  by  a  boundary  condition. 
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APPENDIX  B 


Derivation  of  Distribution  Function 
Boundary  Condition  for  an 
Exponential  Energy  Dependence 

The  purpose  of  this  appendix  is  to  derive  the  expression  for  the 

analytic  source  term,  f  (v),  Eqn  (17).  This  expression  is  valid  if  the 

v 

electron  energy  emission  spectrum  is  assumed  to  be  exponential,  which  is 
a  good  approximation  for  b^ackbody  photon  sources.  The  derivation  given 
here  follows  along  the  reasoning  of  Carron  and  Longmire  (Ref  4),  and 
uses  the  angle  and  direction  definitions  of  Fig.  25. 

For  an  exponential  dependence,  the  electron  energy  distribution, 
dF/dW,  is  related  to  the  X-ray  flux,  <j>  ,  via: 


111 


1 


Since  w  =  total  energy  of  the  electron,  the  normal  component  of 


energy,  w  ,  is, 
2 


w  =  -Lmv2  (B-2) 

z  2  z 


The  transformation 
coordinates  is: 


from  (v  , 
x 


v 

y 


V 

z 


coordinates  to  (  | 


where. 


dv  dv  dv  =  IvIsinG  d9  dd>  dv  =  v2dfl  dv  (B 
x  y  z  y 


dil  =  differential  solid  angle 
=  sing  d9  d^  . 


Therefore, 


dF  d3  F  m  d2  F 

— *  .  - - -  * - 

dv  dv  dv  dv  v  dftdv 
x  y  z 


(B-4) 


where  the  identity,  dw  =  mv  dv,  has  been  used. 

Tow,  the  and  directions  are  integrated  out  of  the 
distribution,  with  the  help  of  Eqn  (B-4),  and 


I  ,  6  ,  ♦  ) 


■3) 


velocity 
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r  z  x  y 


Co  get: 


+“> 

d_F  _  f  f  d3F 

dv  '  /  Idv  dv  dv 
Z  J  J  X  y  z 


dv^dVy 


(B-6) 


since,  dvxdv^=  Vj.dvrd4|. 

The  assumption  is  now  made  that  the  electron  emission  spectrum  can 
be  separated  into  energy  and  angle  parts,  with  a  cos  0  angular 
dependence,  and  exponential  energy  history.  That  is, 


=  MU  cos9  <rwM 

diidw  TTWj 


Yv, t )  e-w/wx 

TTVWj 


(3-7) 


The  1/ir  is  a  normalization  factor  which  is  required  to  give  the  correct 
integration  over  the  hemisphere, 


dF 

dw 


113 


Equation  (R-7)  can  be  substituted  into  Eqn  (B-6)  to  give: 


dF 

dv 


/OO 

v^e 

w 

0 


dv 

r 


(B-8) 


But , 


w 


in(y2 
2  z 


V2) 

r 


And,  dw  =  mv  dv  if  v  is  held  fixed.  As  v  runs  from  0  to  «  . 
r  r  z  r 

from  mv^/2  to  00 .  Thus, 


Using  the  definition  of  the  exponential  integral, 


E  (x) 
1 


dt 


the  differential  flux  is: 


w  goes 
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(B-10) 


dF  _  mYvz  i?(  t )  /_1_  nvj  \ 

dvz  wj  1^2  wt  y 


But,  the  differential  rate 


for  emitting  electrons  at  z=0  is  just: 


dV  =  vzf(0*vz*t)  (3_11) 


That  is  to  say,  the  flux  of  electrons  at  the  boundary  is  equal  to  the 
integrated  flow  rate  across  the  boundary. 


F(t) 


v 


/ 


max 

v?f(0,v 


z 


,t)  dvz 


Comparing  Eqn  (B-ll)  to  Eqn  (B-10), 


f(0,vz,t) 


mY<j>(t)  l _j_  mvj\ 
W1  ^\2  W1  / 


( B— 12) 


Finally,  note  that 


j>(t) 


[sec-1  ] 


C  B— 13) 
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When  Eqns(B-12)  and  (B-13)  are  combined,  with  the  separation  of  Eqn  (46) 
assumed,  the  result  is: 


f 

v 


This  is  the  same  as  Eqn  (17),  when  the  z-subscript  is  dropped. 
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APPENDIX  C 


Irregular  Meshes  and  the  Use  of 
Non-Rectangular  Elements 

Rectangular  elements  were  used  exculusively  throughout  this  disser¬ 
tation  because  of  the  many  simplifications  which  take  place  in  the 
finite  element  equations  for  one-dimensional  SGEMP.  However,  there  are 
some  advantages  to  using  shapes  other  than  rectangles,  particularly  if 
more  than  one  spatial  dimension  is  involved.  Of  course,  if  this  were 
the  case,  a  two-dimensional  problem  would  require  five  independent  var¬ 
iables,  two  space  coordinates,  two  velocity  coordinates,  and  time. 
Then,  there  is  a  large  combination  of  elements  which  could  be  used; 
perhaps  triangles  for  the  spatial  variables,  to  take  advantage  of  their 
ability  to  model  odd-shaped  surfaces  accurately,  and  rectangles  for  the 
velocity  variables,  to  exploit  the  ease  of  integration  over  this  type  of 
element  (see  bel  No  matter  what  combination  is  chosen,  the  fact 
that  rectangles  are  not  exclusively  used  will  make  the  analysis  more 
difficult. 

In  order  to  see  some  of  the  extra  analysis  which  would  be  required, 
this  appendix  will  consider  the  use  of  triangles  in  the  one-dimensional 
problem.  Figure  26  shows  a  triangular  element  in  a  local  area-weighted 
coordinated  system.  The  primary  advantage  of  triangles  over  rectangles 
is  that  one  can  cover  a  region  with  nodes  in  any  desired  location  by 
using  large,  small,  and  elongated  triangles.  A  mesh  showing  this  type 
of  covering  is  given  in  Fig.  27.  The  very  nature  of  the  arbitrary  lo¬ 
cation  of  the  nodes  when  using  triangles  is  the  cause  for  the  more  dif- 
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ficult  analysis. 

The  linear  trial  functions  used  to  approximate  the  unknowns,  E(z,t) 
and  f(z,v,t),  for  the  element  of  Fig.  26  are: 


And, 


T  =  (a  +  b  z  +  c  v)/2 L  (C-l) 

i  i  i  i 

A  =  area  of  triangle  123 


1  xi  n 

1  x2  7Z 

1  x3  ?3 


(C-2) 


2 

V 

z 


i+lVi+2 
i+l~  Vi+2 
i+2  i+1 


i+2Vi+l 

(C-3b) 

(C-3c) 


(C-3a) 


where  the  subscripts,  i,  are  cyclic;  that  is,  i=l 23 1 23 1 . . . 

The  area-weighted  coordinates  for  the  triangle,  (T3  ,  T2  ,  T3  )  are 
obviously  a  dependent  set,  with  T3  +  T2  +  T3=l.  The  transformation  back 
to  the  (z,v)  global  coordinates  is  accomplished  by: 


z  =  Tj  Z|  +  T2  z2  +  T3  z3  (C— 4a) 
v  =  TjVj  +  T2v2  +  T3v3  (C-4b) 

In  this  local  system,  the  linear  trial  functions,  N^,  i= 1 ,2,3  are 
equal  to  the  coordinates  themselves.  That  is,  =  T^ ,  i  *  1,2,3. 
These  trial  functions  are  excellent  candidates  to  approximate  f(z,v,t) 
using  Eqn  (22).  However,  unlike  the  rectangular  trial  functions  of  Eqn 
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(29),  the  triangular  functions  cannot  be  written  as  a  product  of  two 
functions,  each  only  dependent  on  one  of  the  global  coordinates,  like 
Eqns  (32).  This  causes  "coupling"  of  the  z  and  v  directions  in  the  tri¬ 
angle  elements  (or  any  other  non-rectangular  element).  For  example,  in 
Fig.  26,  E(zq  ,t)  is  calculated  by  integrating  the  product  of  v  times 
f(zQ  ,v,t)  along  the  dotted  line,  which  intersects  several  triangles. 
The  upper  and  lower  limits  of  each  of  these  triangles  are,  in  general, 
functions  of  z.  Therefore,  a  large  amount  of  computer  coding  would  be 
required  to  calculate  this  line  integral  on  an  element-by-element  basis. 
Search  routines  would  be  necessary  to  determine  which  triangles  are 
being  crossed,  and,  the  functional  dependence  of  each  upper  and  lower 
limit  must  be  evaluated.  This  probably  would  result  in  numerical  inte¬ 
grations  for  all  element  matrices,  requiring  more  computational  logic. 
For  rectangular  elements,  on  the  other  hand,  the  upper  and  lower  limits 
of  the  element  are  never  functions  of  z.  Furthermore,  the  regular  pat¬ 
tern  of  the  rectangular  elements  allows  the  FEM  to  be  applied  to  Eqn  (6) 
without  having  to  use  the  Method  of  Weighted  Residuals. 

This  is  related  to  the  question  of  separate  meshes  for  the  Vlasov 

Equation  and  Ampere's  Law  raised  in  Chapter  III.  Consider  Ampere's  Law, 

with  the  trial  functions,  M  (z),  J=1,...,N  ,  and  N  (z,v),  i=l,...,N,  for 

J  z  i 

f(z,v,t).  The  N  can  be  taken  as  the  same  functions  used  in  Vlasov's 
Equation  for  this  regular  mesh  scheme.  The  projection  of  the  N  's  on 
the  one-dimensional  z-space  produces  the  same  nodal  placement  as  would 
be  picked  for  a  solely  one-dimensional  problem.  Thus,  Eqn  (6)  can  be 
expressed  as: 

N 

z  N 

^Jmj(z)Ej  -4^  Jv  N^(z,v)  dv  f  i  ( t )  ( C— 5 ) 

J»  1  i*W  all 

v 
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But,  since  any  piecewise  linear  polynomial  will  do  for  the  finite 
element  approximation  for  E,  the  Mj  (z)  can  be  taken  as  the 
one-dimenional  projection  of  the  bilinear  functions,  N  ^(z,v)  onto  z. 
Therefore,  only  the  nodal  values  of  E  are  needed.  They  can  be  evaluated 
at  the  k-th  node  using  the  fact  that  Mj  (zk>  =  6Jk* 

N  f 

Ek  =  fQl=1 JvNi(zk»v)  dv  fi(t)  (C_6) 


Therefore,  in  order  to  calculate  the  electric  field  at  a  node  of  z, 

only  the  integral  on  the  right-hand-side  needs  to  be  evaluated,  not  Eqn 

(27).  However,  if  N  (z,v)  did  not  have  a  z-node  at  z  ,  the  integral  on 

i  k 

the  right-hand-side  of  Eqn  (C-6)  would  be  a  function  of  z  .  Thus,  sim- 

k 

plification  of  this  integral  only  takes  place  is  for  the  rectangular 
elements.  Then,  Simpson's  rule  can  be  applied  directly  to  Eqn  (C-6),  as 
is  done  in  Chapter  III. 
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APPENDIX  D 


The  Assembly  Process 

The  purpose  of  this  appendix  is  to  fully  describe  the  meaning  of 
Eqns  (36)  and  (38)  which  symbolically  indicate  the  relationship  between 
the  local  element  matrix  and  the  global  matrix.  The  following  concepts 
have  been  taken  from  Oden  and  Reddy  (Ref  34,  Chapter  6),  and  applied  to 
the  matrix  assembly  needed  for  the  one-dimensional  SGEMP  equations. 

Let  { x^}  ,  i=l,...,N,  represent  the  set  of  N  global  nodal  points, 
and,  {x^e)}  ,  i=l,...,N^e\  e=l,...,N  ,  represent  the  set  of  local 

nodal  points  for  each  element,  e.  Then,  the  transformation  which  takes 
the  local  nodal  points  into  the  global  nodes  is: 

N(e) 

x.  =  A(e)x(.e)  (D-l) 

i  Z— f  ij  J 

j-1 

where, 

1,  if  node  j  of  element  (e)  coincides 
with  node  i  of  the  global  set; 

0,  otherwise. 

The  [A^e^]  is  a  Boolean  rectangular  (N  x  N^e')  matrix  which  maps  the 

numbering  system  of  the  local  nodes  in  element,  (e),  into  the  global 

numbering  of  the  nodes.  The  set  e  =  I,...,N  ,  is  a  collection 

e 

of  all  mappings,  which  puts  the  nodes  of  all  the  elements  into  their 
proper  location  in  the  global  model. 
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Consider  Fig.  28,  which  shows  a  set  of  three  elements,  one 

rectangle  and  two  triangles.  The  local  and  global  numbering  chosen  for 

the  nodes  are  also  shown.  In  this  case,  N  =  6,  Ng  =  3,  =  3, 

(2)  (3) 

N  =  4,  and  N  =3.  Thus,  for  element  1, 


M 

x2 

x3 

*4 

XS 

-x6. 


— 1 

o 

o 

r  (i>i 

0  1  0 

XI 

0  0  0 

(1) 

0  0  0 

x2 

0  0  0 

(1) 

.0  0  1. 

L  i  J 

x{!) 
x2 
0 
0 
0 


x3 


(1) 


(D-3) 


The  above  6x3  matrix  is  [A^^].  The  other  two  Boolean  matrices, 
(9)  (3) 

[A  "  ]  and  [A  ],  can  be  written  in  a  similar  fashion.  The  set, 

[  A  ^ 1  ^  ] ,  [A^“^  ] ,  [A^]  fully  describes  the  process  of  transforming 
from  the  local  nodes  to  the  global  nodes. 

Using  these  Boolean  matrices,  it  is  now  possible  to  express  in  more 
detail  the  meaning  of  Eqn  (36): 


<A>ij 


titti ye>  -“W  I?5 

(e)  =  l  k-1  fl  J_Y 


(C,n)N^e)(E,n)  d£dn 


E  <D-4) 

(e)-l  k=l  1  =  1 


C  ^ ) 

The  sums  over£  and  k  run  from  1  to  4  because  N  =4.  The  matrices 


1  24 


(e) 

[A  ]  are  dimensioned  N  x  4,  and  there  are  of  them.  The  expression 

for  [A  ]  is  similar: 
z 


(A  )  .  . 
z  ij 


N  .  . 

e  4  4 

EE 

(e)-l  k= 1  A-l 


+1  +1 


^e^*(e)  h(e)  f  f  (e) 

,Aik  Aj*  b  JJ(h  " 


-1  -1 


.  (e). 

+  vc  )  • 

(e)  3N^e^ 

wk  (S>nh^r A  (C,n)  d£dn 


N 

e 

E^ 

(e)  =  l  k=l  1=1 


v  We).  <eE 

>  >  Aik  (aZ 


1,  J 


1 . N 


(D-5) 


The  method  for  transforming  [A^]j  is  governed  by  a  different  rule, 
Eqn  (38).  Like  [A]  and  [Az]  ,  [ A^,] ^  can  be  written: 


+1  +1 


l(Av>j'y 


it  LEW.’  «« 

/  _  \  _  1  1 1  rt  _  1  J  J 


(e)=l  k=l  il  =  l 


-1  -1 


(D-6) 


But  Eqn  (D-6)  can  be  further  reduced  by  expressing  Mj(z)  as  a  function 
of  only  the  local  coordinates.  In  order  to  do  this,  a  one-dimensional 
z-mesh  must  be  defined  which  is  coincident  with  the  z-portion  of  the 
rectangular  (z,v)  grid.  Then,  the  functions,  Lj(^)  and  L  2(£  ),  Eqns 
(31),  can  be  used,  and  M^(z)  can  be  formally  written: 

N  -1  2 


M/z) 


<e2) 

AJn  Ln  *«> 


(D-7) 


(e  )“1  n-1 
z 
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The  (ez)  elements  are  the  linear  segments  between  two  nodes  at  constant 

v  in  the  (z,v)  mesh.  The  Boolean  matrix,  [A  2  ],  dimensioned  2x(N  -1), 

z 

must  be  constructed  to  project  the  two  nodes  for  each  z-segment 

onto  the  global  nodes  of  the  rectangular  mesh.  For  example,  consider 

the  two  rectangles  in  Fig.  29.  Here,  N  =  2,  N  *  6,  N  *  2,  so  that: 
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Now,  Eqn  (D-6)  may  be  written  as. 


N  N  -1 

e  z  4  4  2 


[(A  ) 
v  J 


k  l 


(e)=l(e  )=lk=l  *  =  1  n=l 
z 


where, 


(a(e)) 
Uvn  'k*. 


+1  +1 


‘(e)/*  f L^ez)(C)W<e)(C,n)i^e)(C,n)  dCdn 

J  J  n  k 


-1  -1 


(D-8) 


This  equation  is  the  precise  definition  of  what  is  meant  by  Eqn  (38). 
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APPENDIX  E 


Derivation  of  Element  Matrices 

Equations  (AO)  show  the  end  product  of  performing  the  integrations 
in  Eqns  (37)  and  (39),  which  define  the  A  x  A  element  matrices,  [a], 
[az],  [a^jl,  and  [av?].  Since  all  16  members  of  each  matrix  has  been 
derived  separately,  it  is  not  practical  to  include  them  all  here.  How¬ 
ever,  for  purposes  of  illustration,  one  member  of  each  of  these  matrices 
will  be  derived  to  illustrate  how  the  final  expressions  for  them  are  ob¬ 
tained. 

All  of  the  integrations  over  the  element  rectangle  are  reduced  to 
integrations  over  the  linear  functions,  Li  (x)  and  L2(x),  defined  by  Eqns 
(31).  The  identities  given  in  Table  VI  are  useful  in  these  derivations. 

First,  consider  Eqn  (A7a).  Using  the  definition  given  by  Eqns 
(33d)  and  (32b), 


+1  +1 


ai+2  =  ab 


If 


[MO  +  a3FU)][L2(n)  -  a4F(n)]L2(?)L1(n)  d£dn 


-1  -1 


*  ab 


f+1 

I  [L1(5)L2(C)  -  3ct3L1U)4(£)] 


-1 


f+1 

/  [Li(n)L2(n)  +  3a4L2i(n)L2(ri)]dn 
-1  (E-l) 


where  the  meaning  of  F(£)  and  F(q)  have  been  used,  Eqn  (3A).  Table  VI 
gives  the  values  for  the  integrals  in  Eqn  (E-l),  which  now  can  be 
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TABLE  VI 


Useful  Integrals  of  Lj(x)  and  L2(x) 


+  1 


/L 

-I 


L^(x)  dx  =  1 


+1 


/L 

-1 


L2(x)  dx  -  ~ 


f 

-1 


L  (x)L.(x)  dx  =  ■=■ 
i  3  3 


-1 


L^(x)  dx  = 


I L2 (x)L  (x)  dx  =  4- 

J  i  3  6 

-1 


f 

-1 


1 

L2  (  x  )  L2  (  x  )  dx  =  -L 
i  3  15 


/V>Lj(x)  dx=^ 

-1 


i,3  =  1,2 
i  ^  j 


expressed: 


au2  -  ab((I-  ia3)(I+  -  -|£(2  -  3a3)(2  +  3a,,)  (E-2) 


36' 


Figure  3  shows  that  for  i"4,  j*2,  ?  “~1 ,  C  *+l,  P  *+l ,  and  n  =-l 
Therefore,  Eqn  (40a)  becomes, 


a42  *  ff<-§+  "  a3}  “  ff<2  "  3a4)(2  '  3a3}  (E'3) 
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identical  to  Eqn  (E-2). 

Now,  recall  the  definition  of  (a  ]  given  by  Eqn  (37b).  In  order  to 

z 

calculate  (a  )  ,  the  following  derivative  is  needed: 

z  23 


|^3(C.n)  »  -|l  ^n) 


Thus , 


+  1  +1 


(a  ) 

Z  23 


-|  /  I  ( bn  +  vc)(L2(£)  -  c^FfONLjCn)  +  <*2F(n)]L2(n)  d£dn 
-l-l 

fl 

j/[(2bL2(n)  -  b  +  v^) (L^rjJL^tr))  -  3»2L[(n)L2( i>)) ]dn  • 


/•+1 

•  J  [L2(0  +  3«x1L1(C)L2(5)]  d£  ,  (E-4) 


since n  «■  2L  2  (  n  )  -  1.  Carrying  out  these  integrations  yields: 


(a  ) 

2  23 


.1  3 


1  1. 


«  -  b><i-TZ?2>  +  -f»2>l 


-lo  +  api^-Vc  «-|<.2)]  (E-5) 


With  i»2,  j*3,  which  implies  that  k-1,  p-2,  ?  *+l,  £  ■*+!,  n  “-1,  h  »flt 
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Eqn  (40b)  gives, 


(az23  "  17°  +  °1  )((2  '  3a2  )vc  -  |a2b)]  (E-6) 


which  is  the  same  as  Eqn  (E-5). 

Now,  for  (a  , )  : 

vl  12 


Therefore, 


<%i>12 


+1  +1 


+  «  iFCC) ] [L x(n)  +  F(n)]L2(0  d£dn 


f1  +l 

-j  lLfaW)  ~  3alL21(C)L2(C)J  d?«  j[ L^n)  -  3al(L1(n)L2(n)]  dn 

J-l 


(E-7) 


The  above  equation  yields  upon  integration: 


(avl)12”f(i"^l)(1'-°,,»)  "  ■•^(1  “V<5  -  <E~8> 
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Again,  by  definition,  with  i=l,  j*2,  p=4 ,  k»l,  £/- 1,  £  »+l ,  n /-l , 
n  .=~1,  Eqn  (40c)  is, 


(a  ) 
vl  12 


~  (1-  ~3a-H)[5  -  -|<9  +  3)ct ! ] 


-  ~  o4)(5  -  6o x)  (E-9) 


identical  to  Eqn  (E-8). 

Finally,  consider  the  evaluation  of  (a  )  .  Since, 

v2  33 


$3«-n>  ■IL2(° 


then, 


+1  +1 


(-v2,33-f//"LV5)(L2<» 


-1  "I 


-  a3F(£)][L2(n)  -  0/(11)]  d£dn 


+1  ±1 
jJ[ L32(£)  +  3a3L|(£)L1(£)]  d£*  j [L2(n 


)  +  3a2L1(n)L2(n)]  dn  (E-10) 


Integration  gives: 


(\2>3  -  H  *  -reV0  -  °2>  •  Iff'5  +  3C,3><‘  *  »2> 


(E-ll) 


With,  i-3,  j-3,  k-3,  p-2,  ^*+1,  5^-+!,  ^-+1,  n^-H,  and  £5«  £  j— 1 , 


Eqn  (40d)  yields, 


<*y2>„  -^Tr^r'5  +i<9 ' 3)11 31  -  2o<!  +  ^<5  +  3-3>  <E-'2> 

which  is  correct. 

All  of  the  other  members  of  [a],  [a  ],  [a  ],  and  [a  ]  are 

z  Vi  Vi 

obtained  similarly,  to  give  the  expressions  shown  by  Eqns  (40). 
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